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1 . INTRODUCTION 


(1.1) Background 

The topics of dynamic loading of gear teeth and the 
deflections of gear teeth due to dynamic loads have been 
treated extensively. 

One such work, presented by Cornell and Westervelt [1], 
utilizes an improved version of a model developed by Richardson 
[4]. The model generates the dynamic loads for a meshing gear 
using a cantilever beam with a cam moving along it, simulating 
the engagement and disengagement of the adjacent tooth (see 
Figure 1-1). These dynamic loads are then used in a dynamic 
model of meshing gear teeth where the two gear hubs act as 
rigid inertia and the teeth as variable stiffness springs as 
shown in Figure 1-2. Of significant importance in this 
investigation is the claim made by the authors that the effect 
of variable tooth stiffness is small, changing the dynamic load 
response slightly compared to a system with constant tooth 
stiffness. 

Another dynamic load response algorithm was developed by 
Wang and Cheng [2-3], where they reported that both the dynamic 
load and the induced dynamic response are highly dependent on 
the speed of the moving load. In slow speed regions, the dyna- 
mic load response is composed of a static response which varies 
with the stiffness of the tooth. Superimposed on the static is 
an oscillatory response caused by the excitation of the system 
at the resonant frequency. Wang states that as the speed of 


1 




2 


the moving load increases and the resonant frequency of the 
system is approached, the dynamic load response becomes so 
abrupt that tooth separation occurs. A much smoother response 
is generated when the speed of the moving load is increased and 
becomes out of phase with the system resonance. Here the 
peak response is reduced significantly and actually becomes 
less than that for a static load. Examples of the dynamic load 
variation obtained by Wang and Cheng are included in Figure 
1-3. 

Kasuba [5] presents an algorithm which analyzes spur 
gearing under static and dynamic loading conditions. In his 
analysis, the stiffness of the teeth are determined by solving 
the statically indeterminant problem of multi-pair contacts, 
changes in contact ratio, and meshing gear deflections. in 
general, Kasuba states that to decrease the dynamic load 
response, increased damping and/or contact ratio can be used. 
He also noted that, in a general sense, high contact ratio 
gears have lower dynamic loads than low contact ratio gears. 

Up to now, the discussion of models developed to determine 
the dynamic response of gear teeth has been limited to theore- 
tical cases. Wallace [9] in his investigation, uses finite 
element analysis in conjunction with experimental techniques to 
study the deflection of gear teeth. He subjects a single 
tooth to both Hertzian impact and general dynamic loads, hoping 
to define a procedure for predicting deformation distributions 
due to dynamic loads. The finite element method shows good 
correlation with experimental results obtained using a short 
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cantilever beam subjected to impact loadings at different posi- 
tions . 

Another important contribution to the subject of gear 
dynamics was made by Attia [6]. He studied the effects of 
including the rim when performing a static analysis to deter- 
mine tooth spring constants. He concluded that the stiffness 
of teeth with the rim included is significantly less compared 
to a variable cross-section cantilever beam rigidly fixed to 
the gear body. With the added flexibility, the initial con- 
ditions of two meshing teeth are highly dependent on the 
deflections of the two previously engaging teeth. This fact is 
very significant, as it will definitely affect the type of load 
experienced by the upcoming gear pair. 

Many of the theoretical models used to predict the deflec- 
tions of gear teeth, such as those presented by Cornell [1], 
Wang [2-3] and Kasuba [5], make use of tooth stiffness 
variations obtained from a static deflection analysis. The 
equations of motion are expressed as functions of the load 
position only. 

Nagaya and Uematsu [7] state that because the contact 
point moves along the involute, the dynamic load response 
should be considered as a function of both the position and 
speed of the moving load. In their analysis, they approximate 
the deflections of actual gear teeth due to moving loads by 
modelling the tooth as a tapered Timoshenko beam. They present 
plots of normalized centerline deflections for different moving 
load speeds, and claim that dynamic stiffness variations can be 
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derived from their results. However, as illustrated later, 
this claim turns out to be false. 

In order to make the theoretical developments of models 
used to determine the dynamic response of gear teeth more prac- 
tical, some assumptions are made. One such assumption made by 
the first three authors presented, is that the mass of the gear 
tooth compared to the gear body is small and can be neglected. 
Literature gives no hints to whether this assumption has been 
thoroughly investigated. 

(1.2) Problem Statement 

In this study, two basic problems are investigated. The 
first phase is to determine whether the dynamic response of a 
single spur gear tooth is dependent on the speed of a moving 
load acting on the tooth. 

The second phase is an investigation to determine the 
significance of omitting the inertia of the gear tooth from the 
dynamic deflection model due to the small mass relative to the 
gear body. 

(1.3) Scope of Work 

A model based on involute geometry is developed to automa- 
tically generate a spur gear tooth profile and finite element 
mesh, including the rim, using a minimum of input parameters. 
This model is then used to determine the effects of the speed 
of a moving load on the deflection of a single gear tooth. Two 
constraint configurations are tested; one where only the invo- 
lute profile and fillet regions are allowed to deform, the 
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other with the entire rim included. The results are first 
represented as normalized deflections of the tooth centerline. 
Then the tooth tip deflection time histories are studied for 
the entire load cycle. 

The second phase of the work is to model a meshing gear 
tooth pair using two cantilever beams attached to moveable 
foundation masses. Relative displacements of the foundation 
masses as well as beam deflections are determined for moving 
load speeds bracketing the system resonant frequency. 
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2 . MODEL DEVELOPMENT 


In order to effectively perform a static and dynamic ana- 
lysis of spur gear teeth using finite element techniques, a 
model is needed to automatically generate a tooth profile and 
the accompanying finite element mesh for different size gear 
teeth. Also, the geometry of the tooth should be defined using 
a minimum of parameters corresponding to those most generally 
specified when generating a tooth profile. One such list of 
parameters is: 


Pressure Angle = 8p 

Pitch Radius = RP 

Addendum = AD 

Dedendum = DED 

Circular Pitch = CIRP 

Backlash = BACKL 

Fillet Radius = RF 

Rim Thickness = RTH 


With these parameters, the profile of any spur gear tooth can 
be generated including the rim. 

In the proceeding sections, the equations necessary to 
construct the tooth profile using the preceeding parameters are 
developed, including the implementation of these relationships 
in a profile generation algorithm. The topic of finite element 
mesh generation is also discussed, along with an overview of 
the mesh generation algorithm used to generate the grid. 
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Later in this chapter, a brief discussion of the plane 
strain finite element type used to model the gear is included. 
Also, a general treatment of a linear quadrilateral element is 
used to help develop equations describing the moving loads used 
on the gear teeth. These relations are then implemented in a 
moving load generation algorithm using idealized load time 
history equations for a spur gear tooth. 

(2.1) Profile Generation 

The profile generation sequence is divided into three 
sections; determining relationships, first for the involute, 
then for the fillet, and finally for the rim. 

(2.1.1) Involute Generation 

An involute curve is generated by unwrapping an inexten- 
sible cord from a cylinder. Figure 2-1 illustrates that as the 
cord is unwrapped from the cylinder, point B on the cord traces 
an involute curve AC. The radius of curvature of the involute 
varies continuously, being a zero at point A and increasing 
towards C. At the instant shown, the radius is equal to BE, as 
point E corresponds to the instantaneous center of rotation 
about point B. When generating the involute of a spur gear 
tooth, the cylinder from which the cord is unwrapped corresponds 
to the base circle. This concept is further developed as shown 
in Figure 2-2. Here the local coordinate system X-Y, fixed to 
the hub at the base circle, is used to determine relative coor- 
dinates along the involute. The parameters shown are; the base 
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circle radius RB, the pitch circle radius RP, the roll angle 
0^ and the pressure angle 0p. The pressure angle is defined 
by drawing a line perpendicular to the base circle and passing 
through the point P. This line corresponds to the pressure 
line or line of action of forces between meshing gear teeth. 
The point P being the pitch point. From the triangle OBPO, the 
base circle radius is defined in terms of the pitch radius as: 

RB = RPcos0p (2-1) 

In order to generate discrete points along the involute, 
0i is used in place of 0^ and allowed to vary from zero to a 
maximum, 0 max = ^0i; i =1,2,3, ... ,n, corresponding to the desired 
height of the involute, as shown in Figure 2-3. The maximum 
value of 0 , corresponding to the point B n on the tip of the 
tooth, is found by writing the equation for the triangle 

(RB0max) 2 + R ®2 = R ° 2 
Solving for 0 max gives; 


0 


max 


/ ro 2 -rb 2 \ * 

\ RB 2 / 


(rad) 


( 2 - 2 ) 


where RO is the outer radius defined as the sum of the pitch 
radius and the addendum. Each increment of 0 i produces a point 
on the involute progressively further from the base circle. In 
terms of the local axis system, X-Y, the coordinates of the 
points Bi are determined from the geometry shown in Figure 
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2-4. In simplified form the equations for the X and Y coor- 
dinates are: 


XBj_ = -RB(sin9i - 0j.cos0i) (2-3) 

YBi = RB(cos0i + 0isin0i - 1) (2-4) 


where RB 0 j. is the arc length from the origin of the X-Y system 
to point Ai . 

Next, those equations defining the overall geometry or 
size of the tooth are presented. From Figure 2-5, it can be 
seen that; 

(RB 0P ) 2 + RB 2 = RP 2 


from which; 


•! ' 


(rad) 


(2-5) 


Also from Figure 2-5, it is obvious that; 

>P> 


0 = tan 
P 


■1/RB0JN 

\“rb7 


which can also be written as; 


0p = tan“^0Pj (2-6) 

where 0^ is expressed in radians. During the process of 
calculating actual tooth dimensions, equation (2-6) serves as a 
useful derivational check on qP . With 9 P an ^ 0P t ^ ie an gle 
(p is written as the difference of the two previous angles; 

<J> = 9 P -6p (2-7) 
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Figure 2-5: Spur gear tooth geometry 


And finally, from Figure 2-5, a is found to be; 

a . eP - ta„-l(9?) ♦ 5^™ ( 2 - 8 ) 

where CIRTH is the circular thickness measured on the pitch 
circle, given by; 

CIRTH = CIRP - BACKL (2-9) 

2 

Since one leg of a passes through the tooth center, this angle 
is well suited for transforming the involute coordinates from 
the X-Y axis system to a system whose Y axis passes through the 
center of the tooth, such as Y 1 shown in Figure 2-6. 

When analyzing a gear tooth to determine stresses, deflec- 
tions, etc., it is very advantageous to make full use of the 

axisymmetric properties of the tooth. The involute points 

generated relative to the X-Y axis system are, therefore, 

transformed into another system X'-Y', taking full advantage of 
these properties. 

Using the pitch point on one side of the tooth as a 

reference, as seen in Figure 2-6, a vector f is drawn from the 
pitch point. Bp, to the gear center, ft , which defines the ori- 
gin of the X'-Y' coordinate system. The vector, r', is com- 
posed of two vectors, R and r; 

r ' = R + r (2-10) 

where ; 
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RP 


Figure 2-6: Definition of ax i symmetric 
coordinate system 
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R = d j ' + RBsinai ' 


( 2 - 11 ) 


with ; 


d = RBcosa 


and; 


r = xl + yj (2-12) 

(x and y are the coordinates of point Bp calculated in terms of 
X-Y using equations (2-3) and (2-4)). 

In the new coordinate system, the coordinates of point Bp 
are now defined as; 

X' = RBsin ot + xcos a + ysina (2-13) 

Y' = d - xsina + ycos a (2-14) 

Figure 2-7 illustrates more clearly the elements comprising 
equations (2-13) and (2-14). 

In the profile generation algorithm, included in Appendix 
1, eleven points are calculated along the involute. Equations 
(2-2), (2-3), (2-4), (2-5), (2-8), (2-13), and (2-14) are used 
directly to calculate the point coordinates in the X'-Y' axis 
system. 

(2.1.2) Fillet Generation 

In the present work, two different spur gear tooth 
geometries are considered; low contact ratio gearing (LCRG) 
and high contact ratio gearing ( HCRG ) . By definition, the con- 
tact ratio is the length of the path of contact of mating gears 
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Figure 2-7: Coordinate transformation geometry 
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divided by the base pitch. More practically, it can be thought 
of as the average number of teeth in contact during the meshing 
cycle. A high contact ratio gear is one which has at least two 
teeth in contact at all times. 

One of the main differences between the two forms ( LCRG and 
HCRG) is the fillet transition (see Figure 2-8). For an actual 
low contact ratio gear, the fillet radius is placed tangent to 
the involute and the root circle as shown in Figure 2-8a. The 
amount of overlap of the involute may be different for any 
given tooth design. in the current model, however, the fillet 
radius is calculated to fit tangent to the involute at the base 

ni ml o a fa rman f f f Ko rnr» f 01 ml a c ehr\un in f Ko i f i orl 

A 1. V» »* « WV**4^Si.»a W W\/ fc W W V* | V4 W »»*» *44 W44 W »*4W X* * 4. * 

case (see Figure 2-8b) . 

When designing the high contact ratio gears, the fillet 
region is undercut to provide additional clearance for the 
engaging teeth. Also, the HCRG tooth is generally longer due 
to addendum or other profile modifications, thus the radial 
distance between the base and root circles is also extended as 
shown in Figure (2-8c). 

Given the gear parameters defined for a particular gear, 
the following equation can be used to determine whether the 
gear is a low or high contact ratio gear [10]. 

RR2 + 2RF RR > Rb2 (2-15) 

In equation (2-15) RF is the fillet radius specified for a 
given tooth. If this inequality is satisfied, the tooth is 
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classified as LCRG. This means that the specified fillet radius 
will overlap the involute, and thus must be changed to fit the 
modified form as described in Figure 2-8b. 

Figure 2-9a illustrates the geometry used in developing 
the LCRG fillet radius relations. From Figure 2-9a it is 
obvious that the fillet radius needed to make the transition 
from the base circle to the root circle will have to be larger 
than the radial distance DD as shown. The equation of the 
given triangle is; 


(RF + RR) 2 = (RR + DD) 2 + RF 2 


(2-16) 


By expanding and rearranging terms, equation (2-16) can be 
written for RF as; 


RF = 


2RR DD + DD' 
2RR 


(2-17) 


This then gives the equation of the fillet radius which will 
fit tangent to the involute at the base circle, and tangent to 
the root circle. 

Rewriting equation (2-15) with the inequality reversed 
gives the equation defining a HCRG. For HCRG the transition is 


RR 2 + 2 RF RR < RB 2 (2-18) 

divided between a radial line tangent to the involute and a 
fillet radius from the end of the tangent line to the root 
circle (see Figure 2-8c) . Instead of calculating a new fillet 
radius, as done for LCRG, equation (2-16) is used, along with 
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the specified fillet radius, to calculate the length of the 
radial line DR (see Figure 2-9b). Rewritten in another form, 
equation (2-16) becomes; 

DD 2 + 2RR DD - 2RF RR = 0 (2-19) 


and can be used to determine the radial distance DD spanned by 
the fillet radius. Using the positive root of the quadratic 
equation (2-19) for DD yields; 


DD= 


-2RR+ 


(4RR 2 + 8RF RR) ^ 


( 2 - 20 ) 


DD is then subtracted from the difference between the base and 
root radii to give the length of the radial tangent line. 


DR = (RB - RR) - DD 


( 2 - 21 ) 


When programming the preceeding equations to calculate the 
fillet coordinates, eight equally spaced points are used. For 
LCRG, the arc AOB is divided up into eight equal angles, 6i 
(see Figure 2-10a). Coordinates of successive points are 
calculated by adding 9 i's together for i=l,2...,8 until the 
arc from A to B is generated. The coordinates of Bi in Figure 
2-10b are found from; 


x Bi = X A “ RFcose 
Y Bi ■ Y A " RFsinQ 

For HCRG, the radial distance required for the fillet 
radius, and the radial distance of the tangent line may vary 
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from one gear design to another. To insure equal point 
spacing, integer arithmetic is used to weight the number of 
points between the radial portion and the fillet radius 
according to their respective sizes. This is done to facili- 
tate the finite element mesh generation routine discussed 
later. As in the involute profile generation scheme, the 
points for the fillet are calculated in the X-Y coordinate 
system, and then transformed to the X'-Y' system. 

(2.1.3) Rim Generation 

With the involute and fillet defined, the rim is then 
generated. As stated previously, the fillet radius is placed 
such that it is tangent to the root circle for both LCRG and 
HCRG . From this tangent point, the rim of the gear is added by 
drawing an arc on the root circle. The distance the arc is 
extended on either side of the gear is approximately equal to 
the circular thickness of the gear as shown in Figure 2-11. 

The angle ADA is determined from; 

ADA „ sin-1 ( ™ra i) - sin-iffg) (2-22) 

This angle is then divided into six equal segments and the 
coordinates of points on the rim are calculated using a method 
similar to that shown in Figure 2-10. From the last point on 
the root circle, coordinates for a radial line extending inward 
a distance equal to the specified rim thickness RTH are calcu- 
lated. To complete the tooth profile, coordinate points on the 
inner portion of the rim are calculated using a technique simi- 
lar to that used for the outer rim. 
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(2.2) Finite Element Mesh Generation 

No absolutely correct method has been found to model a 
system using a finite element mesh, even though the topic of 
mesh development has been treated quite extensively. With dif- 
ferent element types and solution techniques, several equally 
valid methods are available for any particular application. 

Indeed, Cook [14] states that although an optimum mesh can 
be determined by requiring that element boundaries follow lines 
of constant strain, this optimum condition only exists for 
one set of loading conditions. As the load changes, so does 
the optimum mesh configuration, and for problems involving 
other than static loading, the difficulties are compounded. 

However, when developing a mesh, simple guidelines can be 
followed which will produce a well enough refined mesh to 
obtain more than satisfactory results. To mention a few; ele- 
ment boundaries should be aligned with structural or geometric 
boundaries and principal load trajectories, element aspect 
ratios should be kept low (less than 7), and when different 
element sizes are used transitions between different size ele- 
ments must be gradual (mesh grading). 

The finite element mesh generation algorithm used for this 
analysis was developed in accordance with the preceeding rules, 
as well as maintaining computational efficiency. Figures 2-12 
and 2-13 show the nodes and elements, respectively, for a low 
contact ratio gear, and Figures 2-14 and 2-15 illustrate the 
high contact ratio geometry and also the varying rim thickness. 

The grid consists of 319 nodes and 276 quadrilateral ele- 
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Figure 2-13: Elements 
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ments. Ten equally divided vertical rows are used to form the 
involute portion of the gear (elements 1-100). Nodes on the 
surface of the involute correspond to the actual coordinate 
points calculated in the profile generation routine. Close to 
the surface of the involute, the element spacing is small pro- 
viding additional stiffness for the application of the load. 
Towards the center of the tooth the element spacing is greater 
where less stiffness is needed. 

The transition from the end of the involute to the root 
circle is accomplished using one of the two techniques 
described in section (2.1.2). For both LCRG and HCRG, eight 
equally spaced rows of elements are used for the transition, 
again using the actual coordinate points calculated in the pro- 
file generation section as surface nodes. When using the HCRG 
transition, with the radial line and fillet radius, the eight 
surface nodes are divided between the two sections keeping 
nodal spacing as even as possible. 

In order to maintain continuity between different gear 
geometry finite element meshes, elements 1 through 204 remain 
the same size relative to the actual tooth sizes. In other 
words, no changes are made in the grid geometry during the 
generation of a particular gear model. The exception to this 
rule is that elements 205 through 276 do vary in size depending 
on the rim thickness. Figures 2-14 and 2-15 show this variation. 

The algorithm containing the equations developed for the 
profile geometry, as well as those relationships used to create 
the finite element mesh is included in Appendix 1. 
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(2.3) Element Description 
(2.3.1) Planar Elements 

Two element types are considered for this analysis; plane 
strain and plane stress. Due to the geometry and loading con- 
ditions of the tooth, it is modelled as a plane elastic problem. 
A plane body is a region of uniform thickness contained within 
two parallel planes. When the thickness of the body is large 
compared to the lateral dimensions, the problem is considered to 
be plane strain. If the thickness is small, it is considered to 
be plane stress. The difference between plane strain and plane 
stress elements is evidenced in the material property matrices. 
For isotropic materials, the material property matrix for the 
case of plane strain is; 


[D] = 


£(1 — H ) 


(1+M)(1-2 M ) 


1 

0 


1 

0 


0 

0 

(1 2 /*)/ 2 ( 1 - n) 


When plane stress exists 


[*] 


E 

1-M 2 


1 

M 

0 


M 0 

1 0 

0 ( 1 -/ 0/2 


where E=30.E6 is the elastic modulus and y=0.3 is Poisson's 
ratio. The matrix multiplication factor is larger for plane 
strain than for plane stress. 


plane strain: 


E( 1-M) 
(1+m)(1-2m) 


4.0385E7 
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plane stress: 


= 3.2967E7 

l -/* 2 


The matrix elements are also larger (except for element 3,3) for 
plane strain. When combined with the strain displacement rela- 
tions to form the stiffness matrix, these differences result in 
an increase in the stiffness for plane strain compared to plane 
stress. 

The thickness of the tooth used in the analysis is 0.25 
inches. Comparatively, the largest and smallest planar dimen- 
sions on the actual tooth (not including the rim) are 0.224 and 
0.081 inches, respectively. Based on the dimensions it is dif- 
ficult to make a judgement on the correct element type for this 
analysis . 

Figure 2-16 shows representative static deflection curves 
for the plane strain and plane stress element types. The addi- 
tional stiffness of the plane strain element is noted. Since the 
difference in deflections between the two element types is small, 
the plane strain element type is chosen for this analysis. 

(2.3.2) General Element Description 

The plane strain element described earlier can be repre- 
sented by a linear quadrilateral element similar to that shown 
in Figure 2-17a. The intersection of the lines which bisect 
the sides of the element form a normalized coordinate system 
, where; 


5 


x 

b 


n 


= I 


a 
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STATIC DEFLECTION CURVES: COMPARISON OF DIFFERENT ELEMENT TYPES 


Figure 2-16: 
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(b) Generalized force application 


Figure 2-17: 
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Between element corners, K and n vary from -1 to +1. The 
displacements within the elements can be written in terms of 
the shape functions for each node, Ni, Nj, Nk and N1 as; 

U = UiNi (x,y) + UjNj (x,y ) + UkNk(x,y) + UlNl(x,y) (2-23) 

where Ui, Uj, etc. define the magnitudes of the displacements. 
If all nodal displacements are zero except for the coefficient 
of Ni(x,y), which is defined as unity, the displacement from 
node i to the other nodes will decrease from unity to zero. 
Using the parameters shown in Figure 2-17b the shape function 
for node i going from i to j is; 


Ni(x,y) = 

Lt 


(2-24) 


where L is the length between nodes i and j in the direction 
of 5 . 


(2.3.3) Moving Loads 

An arbitrary load, P( £, t ), normal to £ is introduced 
whose components are; Px(£,t), Py( £,t) (see Figure 2-17bl The 
effect of the force P(£,t) on node i can be represented by the 
integral of the load times the shape function and thickness 
in the direction of from 0 to L. Component wise; 


Fx^- = At 
Fyi = At 


-U 

/ 


Px( £,t) Ni(x( £) , 


7 


Py(5,t) Ni(x(£), 


y<£))d£ 


(2-25) 

(2-26) 
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where; At is the element thickness (assumed to be unity). 
Inserting the shape function for node i into equations (2-25) 
and (2-26) yields; 



L 



Fxi = 

f Px(£,t) 

<V )d « 

(2-27) 


J 

Or 



Fyi = 

1 Py<5,t) 

( L : 5 )d£ 

Li 

(2-28) 


J 

0 


Equations (2-27) and (2-28) give the load history at node i as 
a function of £(t), resulting from the arbitrary load P(£,t). 

Conversely, the load history at node j is determined by 
considering the shape function obtained when going from node 
j to i with j at zero and i at unity. Here the shape function 
starts at zero and increases to unity as; 


Nj(x,y) = i (2-29) 

Substituting equation (2-29) into equations (2-25) and (2-26) 
yields the force in the x and y directions experienced at node 
j, resulting from P(£,t). 


Fxi = 

Fyi = 


L 

J' Px( t) ( | )d £ 
0 L 

f Py( £,t) ( | )d£ 


(2-30) 

(2-31) 


0 

Equations (2-27), (2-28), (2-30) and (2-31) can be used to 
represent a moving load by introducing the Dirac Delta 
Function. When used in an integral, it translates a given 
function to the origin and gives the value of a function at a 
given time at the origin. The argument of the delta function 
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takes the form of the variation in the position variable. For 
a moving load with constant velocity, the change in position 
is given by the velocity times the time. The arbitrary moving 
load then takes the form; 

P(£,t) = P(t) s( £-Vot ) 

where Vo is the velocity and t the time. Using the delta 
function in the integrand results in all occurences of being 
replaced by Vot. Thus, the four force equations become; 


Fx* = Px(t) 

(L^ 

(2-32) 

Fyi = Py ( t) 

( L— Votj 
L 

(2-33) 

Fx3 = Px(t) 

(5St) 

J-l 

(2-34) 

Fyj « Py(t) 

( ^t ) 

JLf 

(2-35) 


Plotting these equations as a function of time where the magni- 
tude of P(t) is constant, yields to general force histories 
(see Figure 2-18) for a load moving from i to j . 




Figure 2-18: Linear force histories 
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For an actual meshing gear set, the speed of the moving 
load on a single tooth is not constant, but varies linearly 


with time. The time varying speed can be seen to be (see 
Appendix 2 ) ; 

V(t) = RB u) 2 t 

Now, the force equations take a different form with 3 being 
replaced by the displacement resulting from the above velocity; 


where A can replace the quantity RBa>2/2; 

S(t) = At 2 


For the time varying load the force equations then become; 


pi 


F j 


L-At 

P(t) (-X- ) 


At 

P(t) <— l~> 


( 2 - 36 ) 

( 2 - 37 ) 


where the x and y subscripts are assumed. 
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3. DISCUSSION OF NAGAYA ANALYSIS 


Although the problem of theoretically analyzing dynamic 
gear tooth deflections has been treated extensively [1-5] [10], 
models addressing the problem assume that the variation in 
tooth stiffness can be approximated using a static deflection 
analysis. These models assume that the gear hubs act as rigid 
bodies and that the teeth act as variable stiffness springs. 
The stiffness of the teeth varies with the contact position 
along the tooth and is generally arrived at using a static 
deflection analysis such as the one developed by Weber [12]. 
Recently, K. Nagaya and S. Uematsu [7] proposed that since the 
contact point moves along the tooth during the meshing cycle, 
the dynamic load response should be considered as a function of 
both the position and the speed of the moving load. In their 
paper they generate plots of normalized gear tooth centerline 
deflection curves from which they claim the equivalent spring 
constant of gear teeth can be determined. 

(3.1) Approximating A Gear Tooth with a Timoshenko Beam 

In Nagaya' s analysis, the differential equations for a 
tapered Timoshenko beam are written and solved, in the form of 
an eigenvalue problem, from which a modal response analysis is 
used to determine tooth deflections due to moving loads. 
Nagaya assumed a load of constant magnitude moving along the 
beam at a constant velocity from the tip to the base of the 
tooth (see Figure 3-1). 

Using Kara's [8] assumption for the profile of gear 
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L - Ll= 2,25m 
Ho = 2.48m 
Ll/L = .34 

S = (AoL 2 /Io)^= 4.76 


Figure 3-1: Tapered Timoshenko beam 



teeth, Nagaya claimed that the deflections obtained using the 
beam approximation were applicable to any spur gear defined by 
the parameters 


Pressure Angle =20° 

L - LI - 2.25 m 
Ho =2.48 m 
Ll/L = .34 

S = (Aol2/Io) ] ^ =4.76 

where m is the module, L, Ll, Ho are shown in Figure 3-1, and S 
is the slenderness ratio. The module, m, is the pitch diameter 
divided by the number of teeth, measured in inches. When ana- 
lyzing a gear tooth the above parameters are used to describe 
the Timoshenko beam used for the approximation. An example of 
such a comparison is shown in Figure 3-2 where the approxi- 
mating Timoshenko beam is shown superimposed onto the gear 
tooth used in the finite element analysis of Chapter 4. 
Instead of the beam lying tangent to the involute of the actual 
test gear as shown, it should have passed through the tip of 
the involute as illustrated by the inset figure. The inset is 
a correct representation of Kara's assumption for the profile 
of spur gear teeth. This discrepancy, is solely attributable 
to the use of backlash when defining the gear geometry. 
Backlash effectively decreases the width of the tooth. In 
order to better compare the finite element analysis to Nagaya' s 
work, the tooth, when analyzed, is constrained so that only the 
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portion void of interior elements is allowed to deform (see 


Figure 3-2). The foundation and rim are constrained against 
motion. Later on, the deflection of the gear tooth is again 
analyzed with the tooth, foundation, and rim allowed to deform. 

(3.2) Interpretation of Nagaya Results 

When presenting his findings, Nagaya plotted normalized 
tooth centerline deflections versus normalized load position 
for different moving load speeds. Figure 3-3, taken directly 
from reference [7], illustrates these results. The solid cur- 
ves in the figure represent normalized tooth centerline deflec- 
tions, each one for a different normalized velocity, V*. The 
vertical arrows labelled T represent the load position relative 
to the length of the tooth, where X/'L is the ratio of the posi- 
tion of the load on the tooth relative to the total length L. 
The dotted lines are the static curves obtained from the Karas 
analysis. By normalizing these parameters; deflection, load 
position, and velocity, the results then become applicable to 
any size gear tooth. 

Non-dimensional deflections can be represented by; 


W* = AoEW/PL 

where 

Ao = Area of the base of the tooth (in2) 

E = Elastic modulus (Psi) 

W = Actual tooth deflection in 

direction of applied load (in) 

P = Applied load (lbs) 

L = Extended tooth length (in) 
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Dynamic deflection curves for the gear teeth with various 
moving loads 


Figure 3-3 
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The velocity in normalized form is written as; 

V* = V/ jE/p 

where; 

V = Speed of moving load (in/sec) 

E = Elastic modulus (Psi) 
p = Material density (lb/in^) 

In this equation the denominator represents the wave velocity 
in bars. Finally, the position of the moving load is given by; 


T = Vt/ (L-Ll ) 

where; 

V = Speed of the moving load (in/sec) 
t = Elapsed time (sec) 

(L-Ll) = Actual tooth height (in) 


From the plots shown in Figure 3-3, Nagaya claims that the 
deflections of gearteeth, subjected to moving loads, vary with 
the speed of the moving load. That is, for the same values of 
T, the displacements are directly related to the speed of the 
load. He states that for slowly moving loads, the dynamic 
response reduces to the case of a step function impact load for 
small values of T (see T=0.1, V*=0.001 in Figure 3—3). Since 
Figure 3-3 indicates that the dynamic response is dependent on 
the moving load speed (due to effects of inertia forces of the 
mass of the tooth), Nagaya states that the stiffness of the 
tooth must also depend on the moving load speed. he then 
claims that the varying tooth stiffnesses can be determined 
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from these plots. However, as later demonstrated, Nagaya's 
claim that the response, and therefore the stiffness, is depen- 
dent on the speed of the moving load is a false one. 

A major portion of the present work is directed towards 
substantiation of this conclusion. 
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4. FINITE ELEMENT ANALYSIS 


The deflections of single spur gear teeth with moving 
loads acting on them are determined using finite element analy- 
sis. A single gear tooth is used for six different moving load 
cases. First, the same moving load scheme used by Nagaya [7] 
(constant magnitude and speed) is applied to the tooth, which 
is constrained according to the Timoshenko beam approximation. 
Then the load' application on the tip of the tooth is changed 
slightly and the test repeated on the tooth with the same 
constraints. The two preceeding load cases are then applied to 
a tooth allowing the entire model to deform, including the rim. 
Finally, an idealized load function, with variable load magni- 
tude and speed, is applied to the tooth using both constraint 
cases . 


(4.1) Description of Test Gear 

The gear used as the model for this analysis was selected 
at random. The parameters used to define the geometry of the 
gear are; 


0 = Pressure angle 
P 

RP = Pitch radius 
AD = Addendum 
DED = Dedendum 
CIRP = Circular pitch 
BACKL = Backlash 

RF = Fillet radius 


= 20 ° 

= 1.75 (in) 

= 0.125 (in) 

= 0.175 (in) 

= 0.3927 (in) 
= 0.01 (in) 

= 0.05 (in) 
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RTH = Rim thickness = (j . g (in) 
At = Tooth thickness =0.25 (in) 


Figure 4-1 shows the finite element model of the tooth. 
The test gear is a low contact ratio gear (contact ratio = 
1.74) . 

(4.2) Determination of Normalized Plotting Parameters 
and Their Application to the Gear Tooth 

As stated previously, the normalized deflections of the 
gear tooth are calculated using those parameters specified in 
Kara's assumption for the profile of gear teeth. Thus, when 
the deflections are plotted, the only term in the normalized 
deflection equation taken directly from the gear analysis is 
the deflection of the tooth centerline in the direction of the 
applied load which is perpendicular to the centerline of the 
tooth . 

In Chapter 3.1 the equations needed to define the tooth 
profile approximation, according to Karas, are given. The phy- 
sical dimensions, length, area, etc. are defined in terms of 
the module, m. For a standard spur gear the module is defined 
as the pitch diameter per tooth measured in inches, and is 
usually represented by the inverse of the diametral pitch; 

Module = M ■ 1/DP (in) (4-1) 

where the diametral pitch is; 

diametral pitch = DP = 11 = IL—, = 8 (4-2) 

CIRP .3927 
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Figure 4 - 1 : Test gear 
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To further define the test gear, the number of teeth can be 
calculated from; 

number of teeth = N = 2RP*DP = 2 (1.75X8) = 28 (4-3) 

Given either the diametral pitch or the number of teeth, the 
module can easily be obtained. 

m = 0 . 125 ( in ) 

Using the value for the module and the relations of Chapter 
3.1, the dimensions of the approximating Timoshenko beam are 
determined. The height of the corresponding beam becomes; 

( L-Ll ) = 2.25m = 0.28125 (in) (4-4) 

and the extended length; 

L = (~£) = 0.42614 (in) (4-5) 

0*66 

At the base, the beam thickness is; 

Ho = 2.48m = 0.31 (in) 

and thus the area at the base; 

Ao = HoAt = (0.31) (0.25) = 0.775 (in) (4-6) 

Given the above parameters, the non-dimensional deflections can 
be plotted using; 

W* = AoEW/PL (4-7) 

Again, it should be emphasized that although the beam approxi- 
mation (shown in Figure 3-2) does not match the tooth exactly, 
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the dimensions just determined in equations (4-5) and (4-6) are 
still used in the normalized deflection equation when plotting 
results for the actual tooth. 

Nagaya defines the normalized velocity of the moving load 
to be; 

v * - V/Je/’p ( 4 - 8 ) 

and plots four speeds corresponding to V* equal to; 0.01, 
0.005, 0.003, and 0.001. Using equation (4-8), the actual 
velocities, V, are; 



Table 4-1: Actual Velocities 

When the load moves along the involute at a constant speed, the 
values shown in the preceeding table are used directly. 
However, for a meshing gear set, the velocity along the invo- 
lute changes from zero to a maximum velocity, V maXf according 
to equation (A2-6). When the speed of the moving load is 
modelled using this relation, the maximum velocity is defined 
to be the velocity given in Table 4-1. Therefore, the speed 
starts at zero and increases to a maximum speed corresponding 
to those given in the table. By choosing the previously deter- 
mined velocities of Table (4-1) to occur at the base circle 
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where the velocity is maximum, the RPM of the gear for each 
non-dimensional velocity can be determined. From equations 
(A2-6) and (A2-7) the load cycle time for a meshing tooth can 
be shown to be; 


TF = 2S/V max (4-9) 

where S is the distance along the involute from the tip to the 
base circle, and V max is the velocity at the base circle. 
Equation (A2-6) can then be used to calculate the angular velo- 
city of the gear; 

Vmax 

OMEGA = fF*RB (rad/sec) (4-10) 

For each non-dimensional velocity, load cycle times, TF, and 
the RPM's of the test gear are found to be; 


V* 

V m ax (in/sec) 

TF(sec) 

RPM 

0.001 

2025 

.2436E-3 

21470 

0.005 

1012 

. 4875E-3 

10730 

0.003 

607 

. 8127E-3 

6435 

0.001 

203 

. 2436E-2 

2147 


Table 4-2: Variable Velocity Parameters 


To reiterate, the times TF included in Table 4-2 are those for 
which the velocity starts at zero at the tip and increases 
linearly to a maximum value, V max . For a load moving with 
constant velocity, the time for the load to move over the 
involute is simply the distance, S, divided by the velocity. 

Returning now to the non-dimensional parameters plotted by 
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Nagaya, for a constant speed moving load the position of the 
load along the involute is given by; 

T =-¥£_ (4-11) 
i_|— Li J- 

where T varies from 0.0 at the tip, to 1.0 at the root circle. 
A value T=0.8 correponds to a point near the base circle 
radius between nodes 110 and 121 (see Figure 2-12). Equation 
(4-11) is valid only for constant speeds, V. 

(4.3) Description of Dynamic Loading Cases 

The dynamic deflections of single spur gear teeth are 
generated using three loading cases; a constant speed constant 
magnitude load with impact engagement, a constant speed 
constant magnitude load using a finite load engagement rise 
time (for these two loading cases the load is applied normal 
to the tooth centerline), and a load with varying speed and 
magnitude. In this last case the load is applied normal to 
the involute. 

The first of these three loading cases is designed to imi- 
tate exactly the forcing function used by Nagaya. At time 
equal to zero, a load of 1000 lbs, simulating an impact load, 
is applied to the tip of the tooth, and maintained until the 
end of the load cycle (i.e. from the tip to the base circle). 
(Whenever the terms "impact loading" are used, the author is 
describing a step function). To simulate this loading con- 
dition for the finite element analysis, time functions repre- 
senting nodal load histories are calculated for each node on 
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the involute. For the load moving between nodes i and j f the 
force histories are described by equations (2-36) and (2-37); 


F 1 = P(t) (4-12) 

F3 = P(t) (V£) (4-13) 


where; L is the distance between nodes, V is the velocity of 
the moving load, and t the time. With several load value data 
points defined along the involute, the finite element code uses 
these points and linearly interpolates between them to define 
the time functions. The time functions for this loading case 
are shown in Figure 4-2. In Figure 4-2 the node numbers 
correspond to the first eleven nodes on the right involute sur- 
face of the tooth. 

To determine the effect of impact load engagement, another 
test is run using a finite rise time for the load on the first 
node. Instead of the load applied all at once at time equal to 
zero, it starts at zero and gradually increases to the maximum 
value (see Figure 4-3). Here the magnitude of the load is zero 



PER 

2 


t 


Figure 4-3: Finite rise time 
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Figure 4-2: Impact load engagement 






















at t equal to 0.0 and goes up to 3/4 of Pmax by PER/8. The 
rise time is defined as a fixed fraction (PER/8) of the time 
function period. As the speed of the moving load increases, 
the rise time decreases. Therefore, the rise time for V*=0.001 
is ten times greater than for V*=0.01. The time functions for 
this loading case are shown in Figure 4-4. Only the first time 
function is different between Figures 4-2 and 4-4. 

The last loading case tested is one in which the speed 
and the magnitude of the load vary with time. Wallace and 
Seireg [9] give idealized relationships for the magnitude of 
the load on a gear tooth as a function of time and the contact 
ratio. They are; 

P(t) = h Pmax (1- cos ) for: 0.0<t<aTF 

P(t) = Pmax for: aTF < t < (1-a) TF (4-14) 

P(t) = h Pmax (1-cos ) for: (1-a) TF < t <TF 

where; TF is the load cycle time, t is the time along the invo- 
lute, and a is a factor dependent on the contact ratio. A 
value of 0.28 for a is used, corresponding to a contact ratio 
of 1.56. The force history described by equations (4-14) is 
plotted in Figure 4-5. By applying equations (2-36) and (2-37) 
with the load replaced by equations (4-14), the time functions 
generated for this loading case are like those shown in Figure 
4-6. Note that the time function period decreases as the load 
moves down the involute due to increasing speed. 

Equation (4-12) and (4-13), along with equation (4-14) are 
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Figure 4-4 : Gradual, load engagement 












































used in a time function generation algorithm which is included 
in Appendix 1. 

(4.4) Finite Element Test Results 

The results contained in the proceeding sections were 
obtained using the SAP 6 finite element code, implemented on the 
UNIVAC 1100/80 computer facility at Michigan Technological 
University. 

(4.4.1) Comparison of Static Results 

As a preliminary check on the accuracy of the finite ele- 
ment analysis technique applied to gear teeth, static- 
deflections obtained using finite elements are compared to 
those calculated by Nagaya using Kara's assumption for the pro- 
file of gear teeth. Comparisons are made with and without the 
rim included in the analysis. Figure 4-7 shows the plots of 
the normalized centerline deflections obtained using Timoshenko 
beam constraints. The dashed lines are the static deflections 
calculated by Nagaya. From the figure it is apparent that the 
Timoshenko beam (used to produce the dashed lines) is stiff er 
than the tooth. Going back to Figure 3-2, it is seen that the 
beam is considerably larger than the tooth, especially towards 
the base. Thus one would expect the beam to be stiffer. As 
the load is applied closer to the base of the tooth the dif- 
ference between the static deflection curves becomes less 
exaggerated. For T=0.8 the centerline of the tooth actually 
deflects less than the beam. The reasons for this are not 
completely clear. One possible explanation, however, is the 
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fact that as the load is applied closer to the base, the amount 
of local deformation around the point of load application 
increases due to increased nodal spacing. This causes the 
tooth centerline to deform around the local deformation, thus 
decreasing the overall deflection of the gear tooth. (Appendix 
3 includes the actual tooth in the statically deformed con- 
dition, illustrating the increase in local deformation). In 
these figures, the compatibility of the element is not 
violated. The deformation scale factor causes element overlap. 

With the rim included in the analysis, the centerline 
deflections are considerably more severe (see Figure 4-8). 
The curves obtained by Nagaya, represented by the dashed line, 
are exactly those pictured on Figure 4-7. The purpose of this 
set of plots (Figure 4-8) is to emphasize the added flexibility 
afforded by the rim material. (Appendix 3 also contains the 
tooth in the deflected state with the rim included) . 

(4.4.2) Modal Analysis - Determination of Mode Shapes and 
Natural Frequencies 

In Nagaya' s paper, the differential equation for the non- 
dimensional deflection, W* , is derived and then solved numeri- 
cally. The solution to the differential equation (eigenvalue 
problem) includes an infinite number of natural frequencies 
(eigenvalues) and an infinite number of mode shapes 
(eigenvectors). However, he included only the first three 
eigensolutions in the dynamic response analysis. It should 
also be mentioned that the differential formulation is done for 
transverse vibration so only bending modes are included. 
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STATIC DEFLECTION CURVES: COMPARISON OF FEM AND KARAS' TECNIQUES 


Figure 4-8: Rim included 
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For the present analysis the mode shapes and natural fre- 
quencies are determined from the finite element model by 
solving the general equation; 

( [K] - X [M] ) { D} =0 for X = to 2 

for pairs of x and {D>. These results are then used in a modal 
response analysis to determine the response of the tooth to the 
moving loads. For this analysis ten modes are included, with 
both transverse and axial vibration. Appendix 3 includes the 
first few mode shapes and natural frequencies of the tooth for 
both constraint cases. Again, note the difference in flexibi- 
lity between the two models (with and without the rim). 

The first three natural frequencies from Nagaya ' s work 
are compared with those found for the actual tooth. Both 
bending and axial modes are included in the finite element 
analysis, so the first three bending modes from’ this analysis 
are used for comparison (modes 1,3,4) (see Table 4-3). 


NAGAYA FEM 


Mode 

FREQ ( rad/sec ) 

Mode 

FREQ (rad/sec) 

1 

1.092E6 

1 

5.718E5 

2 

3.764E6 

3 

1.436E6 

3 

9.488E6 

4 

2.650E6 


Table 4-3 Comparison of Modal Results 
As expected, the natural frequencies of the beam approximation 
are somewhat higher than those of the tooth (beam constraints). 
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partly due to the additional material toward the base of the 
beam. 

(4.4.3) Dynamic Deflections: Timoshenko Beam Constraints 

(4. 4. 3.1) Impact Loading 

Shown in Figures 4-9a and 4-9b are the normalized cen- 
terline deflections obtained using the impact engagement 
loading case (see Figure 4-2). Results are obtained for load 
positions of T=0.1,0.2, . . . ,0.8. Each solid line represents 
the non-dimensional centerline deflection; 

W* = AoEW/PL 

due to the applied moving load. Remember also that the deflec- 
tions are plotted as a function of position; 



and not as a function of time. So for T=0.1 the solid lines 
show the normalized centerline deflections for the different 
speeds with the load one tenth the distance between the tip and 
the root. Remember also that the time for the load to move 
from T=0 to T=0.1 is different for all velocities, V*. The 
dashed lines correspond to the static deflections obtained for 
the tooth with the load in the position shown. 

Initial examination of these plots suggests that the 
displacements are indeed dependent on the speed of the moving 
load. For slow speeds and low values of T (V*=0.001, T=0.1), 
the deflection is approximately twice the static as Nagaya 
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DYNAMIC DEFLECTION CURVES OF SPUR GEAR TEETH WITH MOVING LOADS 


Figure 4-9a: Timoshenko beam constraints, 
impact load engagement 
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DYNAMIC DEFLECTION CURVES OF SPUR GEAR TEETH WITH MOVING LOADS 


Figure 4-9b: Timoshenko beam constraints, 
impact load engagement 
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claimed. However, these plots do not give an accurate descrip- 
tion of the dynamic deflections. 

A much more representative picture is obtained when the 
actual deflection of the tooth tip is examined over the load 
cycle using small sampling intervals, ( T=0.001 for V*=0.01). 
Figure 4-10 shows the true time history of the tooth tip as the 
load moves from the tip to the base of the tooth. Instead of 
the tooth being in a particular deformed state at load position 
T, dependent completely on the speed, it actually oscillates 
about a datum with constant amplitude and frequency. By 
assuming that the tip oscillates about the static position, the 
amplitude of oscillation is approximately twice the static, and 
is initiated by the impact load at the beginning of the load 
cycle. The only effect the speed of the moving load has is to 
change the number of oscillations per load cycle. Recall that 
the time for the load to move from T=0.1 to T=0.8 is ten times 
greater for V*=0.001 than for V*=0.01. Therefore, approxima- 
tely ten times more oscillations occur for the slower of the 
two speeds . 

From Figures 4-9 and 4-10 we can then conclude that the 
deflections of the tooth do not depend so much on the speed of 
the moving load, but on the tooth position at the instant the 
deflection sample is taken. By roughly lining up the position 
and deflection on the tip deflection curves, the position of 
the tip of the tooth shown in Figures 4-9a and b can easily be 
duplicated . 

The datum about which the tooth oscillates is determined 
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by repeating the previous analysis with the system critically 
damped. The transient caused by the impact load is then 
"filtered" out with only the steady state response remaining. 
Figure 4-11 shows the non-dimensional centerline deflections 
for all four moving load speeds. From the figure it is 
apparent that all centerline deflections lay over the static 
curve. To verify this claim, the tooth tip deflection 
histories are again plotted. In Figure 4-12 the plots show the 
tip following the static curve. 

(4. 4. 3. 2) Finite Engagement Rise Time Loading 

In this test the tooth is subjected to the moving load 
conditions illustrated in Figure 4-4 where the load on the 
first node is gradually applied over a time of PER/8. This 
loading case produces significantly different results compared 
to the impact load test. (Since the normalized centerline 
deflection curves do not accurately represent the dynamic 
deflection phenomenon, they are not included) . Figure 4-13 
gives the tooth tip deflection history for this loading con- 
dition. The tooth still oscillates about the static position 
with constant frequency, but the amplitude varies significantly 
with speed. The reason for this change, from the previous load 
case, is best explained by again considering the load engage- 
ment rise time. 

In Chapter 4.3, Figure 4-3, the rise time is defined as a 
fixed fraction of the time function nodal period (PER/8). For 
V*=0.01, the rise time is ten times less than for V*=0.001. 
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Figure 4-11 
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Figure 4-12 
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It should then be obvious, that as the speed of the load 
increases, the time function for the first node begins to 
approximate the step function impact load of Figure 4-2. 
Examination of the deflection amplitude for V*=0.01 of Figure 
4-13 shows it to be nearly the same as V*=0.01 of Figure 4-10. 

(4. 4. 3. 3) Wallace - Seireg Loading 

For this loading case (time function shown in Figure 4-6) 
both the speed and magnitude of the load vary with time. With 
a contact ratio of 1.56, the load doesn't reach the maximum 
value of 1000 lbs until it is between the second and third 
nodes. Since the magnitude increases smoothly and gradually, 
no abrupt load changes are encountered. 

The tooth tip deflection history curves for this loading 
case are included in Figure 4-14. Due to the nature of the 
speed variation, the deflections are plotted as a function of 
time instead of position as done previously. Here it can be 
seen that the tip of the tooth follows the static deflection 
curve for each of the speeds. This is again due to the slow 
and gradual engagement of the load. 

(4.4.4) Dynamic Deflections: Rim Included 

Including the rim adds flexibility to the system as 
already mentioned. The tooth tip deflection history for the 
impact loading case, shown in Figure 4-15, illustrates this 
fact. As before, the amplitude and frequency of vibration are 
the same for each of the four speeds. However, compared to the 
beam constraint case of Figure 4-10, the amplitude is con- 
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Figure 4-14: Wallace- Seireg loading 
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Figure 4-15 
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siderably larger and the frequency slower. 

When the finite engagement rise time loading case is used 
the same general results are noted as before. That is, as the 
speed of the moving load increases, the rise time approaches 
impact conditions for V*=0.01 (see Figure 4-16). 

Application of the Wallace-Seireg load history equations 
to the rim constraint case produces results which behave 
exactly as before. In Figure 4-17 the tip of the tooth 
deflects in proportion to the magnitude of the applied load. 
The oscillations about the static curve are evident but do not 
contribute significantly to the overall response. 
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Figure 4-16 
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Figure 4-17: Wallace- Seireg loading 
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5. DETERMINING THE EFFECTS OF INERTIA ON THE 
DYNAMIC RESPONSE OF MESHING GEAR TEETH 

In most theoretical models developed to determine the 
dynamic response of meshing gear teeth [1-5], the mass 
(inertia) of the tooth is neglected in the analysis. This 
simplification is based on the claim that the mass of the tooth 
is small compared to the mass of the hub. 

To determine the effect of the inertia of the tooth on the 
dynamic response of a meshing gear pair, a simplified model of 
two cantilever beams attached to foundation masses is used. 

The system analyzed is illustrated in Figure 5-1. 



Two cantilever beams with identical geometric and material pro- 
perties are rigidly fixed to two foundation masses, Ml and M2. 
In this analysis, the masses of the foundations are defined to 
be the same. The two beams are held together by opposing for- 
ces, P, acting on the masses Ml and M2. A massless ball acts 
as the contact point between the beams and moves along the 
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beams at a prescribed speed. As the contact point moves from 
position 1 to position 2 at a velocity V, the change in stiff- 
ness of the beams causes the masses to oscillate in the direc- 
tion of P at the system resonant frequency. To simplify the 
problem, movement only in the direction of P is allowed. 

The dynamic response of the meshing cantilever beams is 
determined for two loading cases; one where the beams are 
assumed massless, the other where the masses of the beams are 
included. The system is analyzed using constant and variable 
speed moving loads of constant magnitude. Both impact and 
smooth load engagement responses are examined by changing the 
initial conditions of the system. These loading cases are ana- 
lyzed using two values for the foundation masses of 1.0 and 
1.0E-4 lbs. 

(5.1) Analysis Using Massless Beams 

Figure 5-2 shows the system used when the beams are 
assumed massless. As the contact point, represented by a 
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massless ball, moves along the beams, the deflection of the 
beams will change due to the variation in stiffness. The 
stiffness of each beam varies with the local coordinate, £ . 
Since the beams are massless, they do not affect the system 
resonant frequency, but follow the oscillations of Ml and M2 
exactly. To determine the dynamic response for the massless 
beam configuration, the differential equations of motion are 
written and solved for this system. 

(5.1.1) Equations of Motion 

The equations of motion for this system are determined by 
first considering each beam-mass configuration as a free body 
(see Figure 5-3). Writing Newton's Second Law as the sum of 



XI *- 

Figure 5-3: Free body diagrams 
the forces acting on each body we have; 

Ml XI » F - P body 1 (5-1) 

M2 X2 = -F + P body 3 (5-2) 

(Here it is assumed that XI and X2 define the positive displa- 
cement direction). Dividing, through by the coefficient of the 
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second derivative term, and subtracting (5-1) from (5-2) 
yields; 


X2 - XI = ♦ P(^+^) 


(5-3) 


To further simplify equation (5-3), the right hand side is com- 
bined and a common denominator is determined. This gives; 


X2 - XI = (P-F) 


(5-4) 


where; 


1 , 1 M1M2 ,1 , 1 x M1+M2 

Ml M2 M1M1 'Ml M2' Ml M2 

Making simple substitutions, equation (5-4) can then be written 
as ; 

MX = P-F (5-5) 

where; 

w _ M1M2 
M1+M2 

and; 


X = X2-X1 


From beam theory, the static deflections of each beam at 
the point of contact of the load, are given by; 


61 = 




3EI 


62 = 


FS, 


3EI 


(5-6) 


as illustrated in Figure 5-4. Equations (5-6) are written in 
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Figure 5-4: Beam deflection configuration 


terms of the beam stiffness as; 

= H ; i2 = iS l5 - 7> 

where; 

KL = S , K = S 

Cl 3 52 

represent the stiffnesses of the beams. 

As the load moves along between the beams, the varying 
stiffness causes Ml and M2 to oscillate at the resonant fre- 
quency of the system. To insure constant contact between the 
beams while the masses are vibrating, the following rela- 
tionship, termed the constraint equation, must be satisfied; 

X2 - XI = 61 + 62 (5-8) 

where 61 and 62 assume orientations as shown in Figure 5-4. 
Substituting equations (5-7) into (5-8) gives the constraint 
equation in terms of the beam stiffnesses as; 


88 



(5-9) 


X2 


XI = 



F 

K2 


Combining the right hand terms in equation (5-9), 

X2 - XI = F(Sjtp.) 


multiplying through by 

v _ K1K2 
K1+K2 

and letting X=X2-X1, yields a familiar form of the constraint 
equation, describing the deflection of a spring; 


F = KX 


(5-10) 


Inserting F from equation (5-10) into equation (5-5) gives the 
equation of motion for the massless beam system; 

MX + KX = P (5-11) 


where; 


X = X2 - XI 
X = X2 - XI 


M1M2 
“ M1+M2 


K 


KIK2 
= K1+K2 


Initially, the system deflection is determined by assuming that 
static equilibrium is satisfied. Thus, at time equal to zero, 
the initial deflection is the applied load divided by the total 
stiffness (the combined deflections of the beams); 
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X(0) = X2 ( 0 ) -XI ( 0 ) = P/K 


With the initial deflections, equation (5-11) is solved 
using a fourth order Runge-Kutta integration algorithm. Since 
the Runge-Kutta algorithm used is designed for systems of first 
order equations, the second order differential equation (5-11) 
must be converted to first order. This is done by defining; 

Y1 = X and Y2 = X 

Substituting these relations into (5-11), we then have; 

MY1 + KY2 = P (5-12) 


where at t=0.0. 


Y1 (0 ) = 0.0 
Y2 (0 ) = P/K 

Given these initial values, the relations; 

Y1 = ( P-KY2 ) /M 
Y2 = Y1 = X 


are used as input for Runge-Kutta and integrated to determine 
the displacement, X=X2-X1 during the load cycle. (See Appendix 
7 for solution algorithm) . 

(5.1.2) Static Analysis 

Before any dynamic analysis is performed, a static 
deflection test is done to provide a reference for dynamic 
deflection comparisons. Solving the relationship; 
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X2 - XI = P/K 


at different load positions produced the plots shown in Figure 
5-5. The three blocks represent the components comprising the 
constraint equation (5-8). Deflections are determined between 
load positions 1 and 2 illustrated in Figure 5-6. The load 
position is measured relative to beam 1 as labelled on the 
abscissa of Figure 5-5. 

(5.1.3) Dynamic Response Results 

The dynamic response of the massless beam system is deter- 
mined for two loading cases; a constant speed 1000 lb load, and 
a variable speed 1000 lb load. For each of these loading 
cases, two sets of initial conditions are considered; those 
defined for static equilibrium in equation (5-11), and another 
set where the initial deflection, X2 - XI, is equal to zero. 
The second of these initial conditions will cause the load P to 
be experienced as an impact load since the beams initially will 
have no deflection and will attempt to return to static 
equilibrium. 

For loads moving with constant speed, speeds of; 1.0, 5.0, 
10.0, 20.0, and 40.0 inches per second are used. This range of 
speeds is chosen in order to bracket the cycle period asso- 
ciated with the fundamental frequency of the system. The fre- 
quency of vibration of the system is constantly changing with 
the position of the load. However, a representative fundamen- 
tal frequency is calculated when the load is halfway through 
the load cycle. At this position the beam stiffnesses are 
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mass 


Figure 5-6: ‘Beam dimensions and 

load cycle definition 
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equal. Then, we simply have a spring mass system whose 
undamped natural frequency is defined as; 


OJ 



_ 1 _ 

2tt 


(cycle/sec) 


With beam dimensions the same as those given in Figure 5-6 and 
the masses Ml and M2 equal to 1.0, the period of oscillation is 
approximately 0.03 seconds. Thus for a constant velocity of 
7.0 in/sec, approximately two oscillations will occur during the 
load cycle. 

When a variable speed moving load is used, the speed ini- 
tially is zero and increases linearly to 1.0, 5.0, 10.0, 20.0, 
and 40.0 in/sec at the end of the cycle. 


(5. 1.3.1) Constant Speed Moving Loads 

For slowly moving loads, relative to the fundamental 
period of the system, one would expect near static deflections. 
This is the case for a speed of 1.0 in/sec as shown in Figure 
5-7. The figure gives the components of the constraint 
equation resulting from an initial deflection equal to the 
ratio of the load and stiffness. The only deviations from the 
static deflection curve are caused by small oscillations at 
the resonant frequency. As the speed of the moving load 
increases, one would expect an increase in the amplitude of 
oscillation as the period of the load cycle approaches the 
resonant period. Increasing to 5.0 in/sec produced larger 
deviations from static due to increased oscillation amplitudes 
(Figure 5-8). This trend continues until the deflection is 
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DYNAMIC RESPONSE OF MESHING MASSLESS CANTILEVER BEAMS 
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maximum for a speed of 20.0 in/sec. (See Figure 5-9 and 5-10). 
In both Figures 5-9 and 5-10, the value of X2-X1 goes negative. 
This means that the beams have separated. Notice also that the 
system at 20 in/sec, does not even begin to approach the ini- 
tial conditions at the end of the cycle, as for the slower 
speeds. A further increase in speed to 40.0 in/sec, shows that 
although the beams themselves deflect appreciably (DELTAl and 
DELTA2), the masses themselves are displaced very little 
because the load cycle is much shorter than the resonant period 
(see Figure 5-11). 

In each case, the system strives toward the static deflec- 
tion position. But due to the inertia of the foundation 
masses, this position may or may not be maintained depending on 
the speed of the moving load. 

As discovered in the finite element analysis (see section 
4.4.3), the type of load engagement significantly affects the 
response of the undamped gear tooth. This is also the case for 
the meshing beams configuration. By changing the initial con- 
ditions of X=X2-X1 to zero, the system reacts to restore itself 
to static equilibrium. Since the system is undamped, a high 
amplitude oscillation is set up due to the rapid movement of 
the foundation masses. This is best illustrated by examining 
Figure 5-12 where the components of the constraint equation 
are plotted for a velocity of 1.0 in/sec. As before, the 
system oscillates about the static deflection position, but 
with very large amplitude. Note also, that separation occurs 
during each oscillation. The same phenomenon occurs for 
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speeds of 5.0 and 10 in/sec (Figure 5-13 and 5-14). However, 
with a further increase in speed, the system does not have suf- 
ficient time to react and the oscillations become less signifi- 
cant (see Figure 5-15 and 5-16). One should not be led to 
believe that the separation occuring for these initial con- 
ditions is caused by the characteristics of the system. It is 
caused, however, by the application of a 1000 lb load to a beam 
which in the physical sense could not support it. However, 
since the system behaves linearly, the characteristics of the 
response are the same for whatever load is applied, and the 
1000 lb load is used simply to exaggerate that response. 

(5. 1.3. 2) Variable Speed Moving Loads 

Introducing a variable speed moving load, versus constant 
speed, has little effect on the general behavior of the meshing 
massless cantilever beam system. The same conclusions can be 
drawn concerning the dynamic resonse trends due to increased 
moving load speed. As a reference, the dynamic response curves 
for both sets of initial conditions and the five different 
speeds, plotted as a function of position, are included in 
Appendix 4 . 

(5.2) Analysis Including Beam Inertia 

Figure 5-17 illustrates the physical system used to deter- 
mine the effects of including the inertia of the beams on the 
dynamic response due to moving loads. This system is exactly 
the same as the one used in Section 5.1 except the mass of the 
beams is included. First, the equations of motion for this 
system are developed. 
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(5.2.1) Equations of Motion 

The equations of motion for this system are somewhat more 
difficult than those developed for the massless beam con- 
figuration. Instead of writing down the differential equation 
directly, a form of Lagange's Equation is used which takes into 
account the added mass of the beams. 

The Lagrange Equation of motion, utilizing Lagrange 
multipliers is; 


d , 3L . 

dt ( 3ql 


_9L_ 

9q± 


+ 


m 

z 

k=l 


X k a ki 


for: i = 1, 2, . . . , n 


(5-13) 


where; 
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L = T - U = Lagrangian 
T = System kinetic energy 
U = System potential energy 
= Generalized coordinates 
Qi = Sum of non-potential forces 
Xfc = Lagrange Multipliers 


a ki 


- = first partial 


derivative of the constraint equation 
with respect to the generalized coordinates 


The last term on the right hand side of equation (5-13) repre- 
sents the constraint forces which makes it possible to regard 
all generalized coordinates, qi, as independent. 

Four generalized coordinates are used to describe the 
system position, velocity, and acceleration. They are; 


XI - describing Ml 
X2 - describing M2 
ql - describing beam 2 relative to 
foundation mass Ml 
q2 - describing beam 2 relative to 
foundation mass M2 


From the definition of generalized coordinates, these four 
quantities are assumed independent of each other, and when used 
together they describe the state of the dynamical system at any 
time. 
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The static deflection of a loaded cantilever beam is given 
by the applied load divided by the stiffness. When considering 
the dynamic deflection of a beam which has mass, this rela- 
tionship is not valid. Instead, the superposition of the 
natural vibrational modes is used to describe the shape of the 
deformation, and when multipled by the generalized coordinate, 
q^, the actual deflection is obtained. 

The mode shapes are an infinite set of eigenvectors 
derived from the differential equation of the cantilever beam. 
To exactly duplicate the deflection obtained from elasticity 
theory, an infinite number of modes must be used. In practice, 
however, only the first few contribute significantly to the 
overall deflection, such as those illustrated in Figure 5-18. 
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Figure 5-18: Natural Modes 
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In this analysis, only the first mode shape is used, the 
natural mode described by the equation; 

<p(£) = C[ (sin3L-sinh$L) ( sin 8£-sinhB£ ) 

+ (cos gL+coshBL) ( cos B£-cosh 30 ] (5-14) 

where the parameters are defined in Figure 5-18. 

We are now ready to develop the equations of motion for 
the mass-beam system of Figure 5-17. The kinetic energy of the 
system is determined by considering the velocity of the foun- 
dations Ml and M2, and the vibration of the beams with respect 
to the foundation masses. By defining the beam deflections to 
be positive as shown in Figure .5-19, the kinetic energy of the 



configuration 


system can be written as; 


T = |m1 XI 2 + ^M2 X2 2 + j p f (XI + (f>i q 1 ) 2 d^ 1 
1 X*2 . . 

+ j p f (X2 - <p 2 q 2 ) 2d Si 

J n 


no 



The potential energy consists only of the strain energy due to 
the beam deflections; 


U 




(5-16) 


0 0 

Where <J> " is the second derivative of the mode shape with 
respect to the length variable, E, . In equations (5-15) and 
(5-16) the mass per unit length, p , the elastic modulus, E, 
and the moment of inertia, I, are constant over the length of 
the beams and can be left outside the integrals. Substracting 
(5-16) from (5-15) to form the Lagrangian we have; 


L 


T - U = jMl XI 2 + |m2 X2 2 





1 f L2 2 

2 EI J 0 ( * 2 "^2> d? 2 


(5-17) 


For the meshing cantilever beam system there are four 
unknown parameters, XI, X2 , ql, q2 and their derivatives, each 
dependent on time, describing the dynamical system. At any 
time, t, the position of the system can be described by the 
constraint equation; 

X2-X1 - 4> 1 q 1 — <}> 2 q 2 = f (Xl,X2,ql,q2) = 0 (5-18) 

If the constraint equation is not satisfied, separation occurs 
between the beams at the point of contact. 

A Lagrange Equation of motion is written for each degree 


in 



of freedom of the system corresponding to the four generalized 
coordinates. This gives four equations relating the genera- 
lized coordinates and X . Including the constraint equation or 
its derivative gives a fifth equation necessary to determine 

In terms of XI, the equation of motion, term by term, 
becomes ; 

= M1X1 + P f (XI + (J) 1 q 1 ) 2 d5 1 
J 0 



3L , 

— i — ) 

3X1 


Ml XI + p 


fd i ♦ 
^0 


q 


i )d ?i 


(5-19a) 


and; 


3L 

3X1 


0.0 


( 5-19 b) 


On the right hand side, differentiation of the constraint 
equation with respect to Xl yields; 


3 f 
3X1 


- 1.0 


(5-190 


Combining equations (5-19a) through (5-190 along with the non- 
constraint force, P, acting on Ml, in the form of equation 
(5-13) we have; 


The term 
defined 


(Ml + p 




as MB1 . 


1 i 2 

d^ 1 )Xl + p r (<() 1 q-^d^ = P-A (5-20) 

^0 


is simply the total mass of beam 1 which is 
Similarly, with X2 as the generalized 


coordinate; 
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P + X 


(5-21) 


f>u2 

(M2 + MB2 ) X2 - p j <^ 2 q 2 d^ 2 = 

Performing the required differentiation of the Lagrangian and 
constraint equation with respect to ql yields; 


3 L 

3qi 


_d_ 

dt 




3 L 

3qi 



XI + (jjj 2 
XI + c )^ 2 

I 2 


^l )d? l 

<i 1 )d £ 1 


1 


( 5-22a ) 
(5-22b) 


3f _ A 

3q x 9 1 (5-22c) 

Combining equations (5-22a) through (5-22c), and noting that no 
external (non-constraint) forces act on the beams, the equation 
of motion derived with respect to ql becomes; 


r L1 

r L1 ? i 

rLl 



p Jo ((j^XDd^+p J 

[ 0 ( *1 

UV 2 
0 1 

q l )d? l = 

(5-23) 

Likewise for q2; 





r L2 .. 

rl>2 2 

fL2 2 



-p J (<D 2 x 2 )d5 2 + P 

^0 J 

o <*2 q 2 )d? 2 +EI j 

'o z 

q 2 )d? 2 = “*2 X 

(5-24) 


The fifth equation of motion is determined by differentiating 
the constraint equation twice with respect to time. This gives 
an equation relating the second derivatives of the generalized 
coordinates with respect to time. Differentiating equation 
(5-18) once yields; 
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After a second differentiation and rearrangement; 

X2-X1- qi~ <J>2 52 = ^1 <31+ 4>2 <32 +2 ( ^1 <31+ $2 32 > (5-25) 

Equations (5-20), (5-21), (5-23), (5-24) and (5-25) are the 

equations of motion describing the dynamical system of Figure 
5-16. 

These equations can be greatly simplified by making the 
following substitutions; 


PHI1 = ; PHI2 = <j> 2 


r L1 

IPHIl = p 4> 1 d5 1 


r L2 

IPHI2 = p 4> 2 dC 2 
^0 


P 



1 



1 


('Ll pL2 

I2PHI1 = El J 4>^ 2 d^ 1 ; I2PHI2 = El 4>£ 2 d£ 2 

F(q 1 ,q 2 ,<l> 1 ,<l>2 ) = $1 q l + *2 q 2 + 2( ^1 q l + ^2 q 2 ) 

(see Appendix 5 for developments concerning these 
simplifications). The equations can then be written in matrix 
from in terms of the second derivatives and as; 


[A]{X> = { B } 


(5-26) 


where; 



(M1+MB1) 


0 

IPHI1 

0 

1 

0 


M2+MB2 

0 -IPHI2 

-1 

IPHI1 


0 

1 

0 

PH 11 

0 


-IPHI2 

0 

1 

PHI 2 

-1 


1 

-PH 11 

-PH 12 

0 

\ 

«" > 
X2 i 

! 


r p 



XI 



p 


(X) = < 

ql 

} and 

{B} = < 

-I2PHI*ql 


q2 



) -I2PHI2*q2 

1 

' x I 

S. J 

1 

1 

I F ( ql f q2 

•s. 

,<I>1,*2) 


The only differences between equations (5-26) and the equation 
of motion for the massless beam system (equation (5-11)) are 
the inertial terms. These are; X1*IPHI1, X2*IPHI2, and the 
terms involving ql and q2 . By eliminating the inertial terms 
from matrix [A], the resulting equations describe the massless 
beam system. Appendix 6 includes the analysis and results from 
this test. 

The initial conditions governing the system described by 
equations (5-26) are chosen such that the massless beam system 
is emulated. When chosing the initial conditions, there are 
two sets of parameters which must be considered. The first set 
consists of the four position variables XI, X2 , ql and q2 . In 
section (5.1.1) it is stated that the beams assume a static 
orientation at the beginning of the load cycle. Thus, by 
chosing ql and q2 as; 


ql ( 0 ) 


P*PHIl 
" I2PHI1 


q 2 ( 0 ) 


P*PHI2 

I2PHI2 


(5-27) 
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and XI (0) equal to zero as before, the deflection X2(0) is 
determined from the constraint equation as; 

X2 ( 0 ) = ql ( 0 ) PHI1 + ql ( 0 ) PHI2 (5-28) 

This same reasoning can be applied to the first derivatives of 
the generalized coordinates, XI, X2 , ql and q2 . Previously, 
the initial velocities of the foundation masses were determined 
such that X2-X1=0 at time equal zero. For this case, the same 
effect is obtained by forcing; 

XI (0) = 0.0 

X2 ( 0 ) = 0.0 (5-29) 

ql ( 0 ) = 0.0 

Then, using the first derivative of the constraint equation the 
value of q2(0) is determined from; 

X2 - XI = q^PHIl + q 2 * D1 PH 1 1 + q 2 *PHI2 + q2*DlPHI2 

where at time equal to zero, all terms are zero except for 
those in equation (5-30); 

q 2 (0) = ( _ q 1 *DlPHIl - q 2 * D1PHI2 ) /PHI2 (5-30) 

where DlPHIl, D1PHI2 and PHI2 are evaluated at the beginning of 
the load cycle and ql and q2 are taken from equations (5-27). 

With the initial conditions described by equations (5-27) 
through (5-30), the dynamic response of the system described in 
equations (5-26) is determined. 

For each time step during the load cycle, the system of 
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equations (5-26) is inverted to determine the column matrix 
{ X} = { Xl X2 ql q2 X} T in terms of the current values of { 1} 
on the right hand side. At the beginning of the load cycle {X} 
is determined for those initial conditions set forth pre- 
viously, along with the right hand side (B). The values of { X} 
are then used with the first derivatives of the generalized 
coordinates (initially equations (5-29) and (5-30)) as input 
for a Runge-Kutta integration routine which integrates for the 
desired solution parameters XI, X2 , ql and q2 . It also deter- 
mines the new XI, X2, ql and q2 which are in turn used for the 
next integration step. The algorithm containing the 

Runge-Kutta (RKGS) and matrix inversion (SIMQ) subroutines 
along with the parameters describing equations (5-26) is 
included in Appendix 7. Therein lies a detailed description of 
the programming steps for the solution of equations (5-26) 
using the aforementioned initial conditions. 

(5.2.2) Dynamic Response Results: Foundation Mass = 1.0 lbs 

As stated in section (5.2.1), only the first material mode 
of vibration is included when determining the dynamic response 
of the meshing cantilever beams where the inertia of the beams 
is included. One would then think that the system whose dyna- 
mic response is composed of a single mode shape would be 
somewhat stiffer than the same configuration where the beam 
deflection is determined from elasticity theory. However, com- 
paring Figures 5-7 to 5-11, with A6-1 to A6-5, it is seen that 
the difference in deflections is negligible. 
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In this analysis, with the beam inertia included, the 
dynamic response is compared to the massless beam problem. At 
slow moving, constant speed loads, the dynamic response con- 
sists of quasi-static deflection with the moving load response 
superimposed over it (see Figure 5-20). (Results for moving 
load speed of 1.0 in/sec were not obtainable due to lack of 
convergence in the Runge-Kutta integration routine). 
Continuing to increase the moving load speed causes an increase 
in the response as the resonant frequency is further excited as 
shown in Figure 5-21 and 5-22. Increasing the moving load 
speed above the system resonant frequency produces a much 
smoother response as illustrated in Figure 5-23. Now, examina- 
tion of Figures 5-20 through 5-23 compared to Figures 5-7 
through 5-11 shows that for values of the foundation masses of 
1.0 lbs and moving load speeds of 5.0, 10.0, 20.0 and 40.0 
in/sec, the responses of the two cantilever beam systems are 
essentially the same. The same conclusions are drawn when exa- 
mining Figures 5-24 through 5-27 as compared to Figures 5-13 
through 5-16. Here, the initial conditions are changed such 
that the beams experience an impact load at the beginning of 
the load cycle. From the figures it can be seen that the 
responses of the two systems are again, very much the same. 

As done for the massless beam configuration, the speed of 
the moving loads is allowed to vary linearly from zero to a 
maximum value during the load cycle. This analysis serves as a 
useful check since the initial conditions are easily defined 
due to a stationary load. The results from this analysis are 
contained in Appendix 4 . 
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DYNAMIC RESPONSE OF MESHING CANTILEVER BEAMS, WITH INERTIA (VEL^CONST) 


Figure 5-20 
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Figure 5-21 
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Figure 5-23 
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Figure 5-24 
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Figure 5-25 
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Figure 5-27 
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(5.2.3) Dynamic Response Results: 

Foundation Mass = 1.0 E-4 lbs. 

In order to substantiate the conclusion that the dynamic 
response is not dependent on the mass of the beam, both the 
massless beam and inertial beam configurations are analyzed 
again with a decrease in the difference between the values of 
the beam and foundation masses. Of course the mass of the beam 
is not included in the massless beam analysis. Previously the 
mass of the foundation was 1.0 lbs and the total mass of the 
beam was approximately 3.6E-6 lbs. In this analysis, the 
foundation mass is decreased to 1.0E-4 lbs while the mass of 
the beam remains the same. Since the overall mass of the 

system is greatly decreased, the corresponding resonant fre- 
quency is much larger. Using an approximation of the resonant 
frequency using a single mass and the load at the midpoint of 
the beam yields an approximate resonant frequency of; 


0 ) 


=71 => 7f?i^ ( ^r ) - 3500 c Y cie / sec - 


and a cycle period of approximately 2.875E-4 sec. Therefore, 
in order to provide an adequate range of moving load speeds 
such that the system resonance is bracketed, speeds of 100., 
500., 1000., 1500. and 2000. in/sec are used. (Approximately 
one oscillation cycle occurs at 1700 in/sec) . To facilitate 
comparisons of respective systems, the results are presented in 
a parallel fashion. Figures 5-28 and 5-29 show the respective 
responses for a moving load speed of 500 in/sec. As seen from 
the plots, the responses are virtually the same. The results 


are similar for the other moving load speeds (see Figure 5-30 
through 5-35). 
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Figure 5-33 
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Figure 5-34 
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6 . CONCLUSIONS 


The objectives of this analysis were twofold. First, the 
dynamic response of a single spur gear tooth subject to moving 
loads was studied to determine the effect of the speed of move- 
ment of the load. A spur gear tooth, modelled using finite 
element techniques, was used in the analysis. From this analy- 
sis it was found that the dynamic response of a single gear 
tooth is not dependent on the speed of the moving load, but 
rather on the type of load engagement experienced at the 
beginning of the load cycle. Including the rim in the analysis 
added flexibility to the system but did not change the overall 
response . 

The second objective was to determine whether or not the 
mass (inertial forces) of the tooth can be neglected when small 
compared to the mass of the gear hub. A simplified analysis 
using meshing cantilever beams was used. For the range of 
speeds tested, it was found that the inertia forces of the 
tooth are small and therefore, the mass of the tooth can be 
neglected when determining the dynamic response of a meshing 
gear system. 
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ooooooooooooooooooooonooon 


*************X*XXXXXXXXXXXXXXXXXXXXXXKA««... A 


* * 

* SPUR GEARTOOTH PROFILE GENERATOR * 

* * 

* WRITTEN BY LYLE W. SHUEY * 

* * 


* THIS PROGRAM GENERATES A FINITE ELEMENT * 

* MESH FOR EITHER A HIGH CONTACT RATIO * 

* GEAR OR A LOW CONTACT RATIO GEAR, IN- * 

* CLUDING THE FOUNDATION AND RIM IN THE * 

* MODEL. * 

* A MINIMUN OF INPUT PARAMETERS ARE NEEDED* 

* TO DESCRIBE THE GEAR GEOMETRY. * 


* 

THEY ARE; 


* 

* 

PRESSURE ANGLE 

: PANG 

* 

* 

PITCH RADIUS 

:RP 

* 

* 

ADDENDUM 

: AD 

* 

* 

DEDENDUM 

: DED 

* 

* 

CIRCULAR PITCH 

: CIRP 

* 

* 

BACKLASH 

: BACKL 

* 

* 

FILLET RADIUS 

:RF 

* 

* 

RIM THICKNESS 

:RTH 

* 

* 



* 


******************************************* 


C****** THE PROGRAM READS THE INPUT VARIABLES 
C****** USING THE PROCEEDING STATEMENTS. 

C 

DIMENSION DUMY (20) ,YPRF(20) ,XPRF(20) 
DIMENSION XX (40) ,YY(40) ,XN(400) ,YN(400) 
DIMENSION XANG (9,9) , XADA ( 9 ) , XINC ( 9 ) 

PI=3. 141592654 
READ ( 5 , * ) PANG 

WRITE ( 6 , * ) • PRESSURE ANGLE : • , PANG 

PANG=PANG * PI/ 1 8 0 . 

READ ( 5 , * ) RP 

WRITE (6,*) 'PITCH RADIUS: ',RP 

READ ( 5 , * ) AD 

WRITE (6,*) 'ADDENDUM: ' ,AD 

READ ( 5 , * ) DED 

WRITE (6,*) 'DEDENDUM: ' , DED 

READ (5,*) CIRP 

WRITE (6,*) 'CIRCULAR PITCH: ',CIRP 

READ ( 5 , * ) RF 

WRITE ( 6 , * ) 'FILLET RADIUS: ' ,RF 

READ (5,*) BACKL 

WRITE ( 6 , * ) ' BACKLASH : • , BACKL 

READ ( 5 , * ) RTH 

WRITE (6,*) 'RIM THICKNESS: ' ,RTH 

RB=RP*COS (PANG) 

THP= ( (RP*RP-RB*RB) / (RB*RB) ) ** . 5 
CIRTH=CIRP/2 . -BACKL 

WRITE ( 6 , * ) ' CIRTH : • , CIRTH 

AL=THP-ATAN (THP) +CIRTH/ (RP*2 . ) 

DP=PI/CIRP 

WRITE (6,*) 'DIAMETRAL PITCH: ' , DP 

RR=RP-DED 

RO=RP+AD 
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onow ooooooooo 


H-4 / 

******************************************* 

* THIS SECTION CALCULATES POINTS ON THE * 

* INVOLUTE SURFACE, FOUNDATION, AND RIM * 

* WHICH WILL LATER BE USED TO GENERATE * 

* NODAL COORDINATES. (BOTH HCRG AND LCRG) * 
******************************************* 

****** GENERATING POINTS ON THE INVOLUTE PROFILE 

DO 5 1=1,11 
IA=1+ (1-1) 

HH= (RO-RB) /10 . 

H=RO-HH*(I-l) 

ANG= ( (H*H-RB*RB) / (RB*RB) ) **.5 
X=-RB* (SIN (ANG) -ANG*COS (ANG) ) 

Y=RB* (COS (ANG) +ANG*SIN (ANG) -1) 

ANGF=ATAN (SIN (ANG) /COS (ANG) ) -AL 
ANGF=ANGF*180 ./PI 

XP=RB*SIN ( AL) +X*COS (AL) +Y*SIN (AL) 

YP=RB*COS (AL) -X*SIN (AL) +Y*COS (AL) 

XX (IA) =XP 
YY(IA)=YP 

WRITE ( 6 , * ) 'XP: • ,XP, 'YP: ' ,YP, 'THETA: ' , ANGF 
CONTINUE 

****** CHECK IF GEAR IS LCRG OR HCRG 

CHECK=RR*RR+2 . *RR*RF 
ITAN»2 

IF (CHECK. LT.RB*RB) ITAN=1 
IF(ITAN.NE.l) GO TO 15 
C 

C****** STRAIGHT LINE TANGENT TO INVOLUTE (HCRG) 

C 

DD= ( -2 . *RR+ ( 4 . *RR*RR+4 . *2 . *RF*RR) ** . 5) /2 . 

D=RB- (RR+DD) 

DDD=RB-RR 
XTAN= ( D/DDD) 

XFILL=(DD/DDD) 

LTAN=INT (XTAN*8 . ) 

LFILL=INT (XFILL*8 . ) +1 
WRITE (6, *) 'DD: ' , DD 

WRITE (6 , *) 'D: ' , D 

WRITE (6,*) ' DDD: ' , DDL 

WRITE (6,*) ' XTAN : ' , XTAN 

WRITE (6,*) ' XFILL: ',XFILL 

WRITE (6,*) ' LFILL: ' , LFILL 

WRITE (6,*) ' LTAN : ' , LTAN 

DRI=D/LTAN 
DR=0 . 0 

DO 10 1=1, LTAN 
IA=12+(I-1) 

XP=XP- ( ( (RB-DR) *SIN (AL) )-( (RB-DRI) *SIN(AL) ) ) 
YP=YP- ( ( (RB-DR) *COS (AL) ) - ( (RB-DRI) *COS (AL) ) ) 
DR=DRI 

DRI=DRI+D/LTAN 
XX ( IA) =XP 
YY (IA) =YP 
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10 CONTINUE 
C 

C********* FILLET RADIUS SECTION (HCRG) 

C 

BETA=ATAN( (RR+RF) /RF) 

BETA=BETA/ LFI LL 
DO 11 1=1 , LFILL 

IA= ( 12 +LTAN ) + ( I - 1 ) 

X=RF-RF * COS (BETA* I) 

Y=-RF*SIN(BETA*I) 

XP= (RB-D) *SIN (AL) +X*COS ( AL) +Y*SIN ( AL) 

YP= (RB-D) *COS (AL) -X*SIN (AL) +Y*COS (AL) 

XX (IA) =XP 
YY(IA)=YP 

11 CONTINUE 
GO TO 21 

C 

C ****** GENERATING POINTS FOR FILLET RADIUS ( LCRG ) 
C 

15 DD=RB-RR 

RF= ( (DD*DD) + (2 . *DD*RR) ) / (2 . *RR) 

WRITE (6, *) 'NEW FILLET RADIUS ' ,RF 
BETA=ATAN ( (RR+DD) /RF) 

WRITE (6,*) 'BETA: ' , BETA*180 ./PI 
BETA=BETA/8 . 0 
DO 20 1=1,8 
IA=12+ ( 1-1) 

X=RF-RF*COS (BETA* I) 

Y=-RF*SIN ( BETA* I) 

XP=RB*SIN (AL) +X*COS (AL) +Y*SIN ( AL) 

YP=RB*COS (AL) — X*SIN (AL) +Y*COS ( AL) 

XX (IA) =XP 
YY (IA) =YP 

WRITE (6,*) 'XP: ' ,XP, ' YP: ' , YP 

20 CONTINUE 
C 

C****** GENERATING POINTS ON OUTER RIM SURFACE 
C 

21 ADA=ASIN ( (XP+CIRTH) /RR) -ASIN (XP/RR) 

DADA=ADA/6 . 

ADA1=ASIN (XP/RR) 

WRITE (6,*) 'ADA: ' ,ADA*180./PI 
DO 30 1=1,4 
IA=20+ (1-1) 

DUMY (I) =10 . 

XP=XP+RR* (SIN (ADA1+DADA) -SIN (ADA1 ) ) 
YP=YP-RR* (COS (ADA1) -COS (ADA1+DADA) ) 

XX ( IA) =XP 
YY ( IA) =YP 

WRITE (6 , *) 'XP: ' ,XP, 'YP: ' ,YP 
DADA=ADA/4 . 

ADA1=ADA1+DADA 
3 0 CONTINUE 
C 

C****** GENERATING POINTS ON RADIAL PORTION OF RIM 
C 

DELTA=ASIN (XP/RR) 

RIM=RR-RTH 

DRTH=RTH/4. 
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DX=RR*SIN ( DELTA) -RIM*SIN ( DELTA) 

DDX=DX/4 . 0 
DO 40 1=1,4 
IA=24+ ( I — 1.) 

XPRF ( I ) =XP-DDX 
YPRF(I) =YP-DRTH 
XX ( IA) =XPRF ( I ) 

YY ( IA) =YPRF ( I ) 

DRTH=DRTH+RTH/4 . 

DDX=DDX+DX/4 . 0 

WRITE ( 6 , *) 'XP: ' , XPRF ( I ) , 'YP: • , YPRF(I) 
40 CONTINUE 

XP=XPRF ( 4 ) 

YP=YPRF ( 4 ) 

C 

C****** GENERATING POINTS ON INNER RIM SURFACE 
C 

DELTA=ATAN (XP/YP) 

XLEG= ( (XP*XP) + ( YP*YP) ) ** . 5 
DELTAl=DELTA/4 . 0 
DO 50 1=1,4 
IA=28+ (I — 1 ) 

DELTA=DELTA-DELTA1 
XPRF ( I ) =XLEG*SIN ( DELTA) 

YPRF ( I ) =XLEG*COS ( DELTA) 

XX ( IA) =XPRF ( I ) 

YY (IA) =YPRF (I) 

WRITE ( 6 , *) 'XP: ' , XPRF (I) , 'YP: ' ,YPRF(I) 
50 CONTINUE 

C 

C ***************************************** 

C * THIS SECTION USES COORDINATES FROM * 

C * THE PREVIOUS SECTIONS TO CALCULATE * 

C * NODAL NUMBERS AND COORDINATES . THERE * 

C * ARE 319 NODES USED IN THIS MODEL. * 

Q ***************************************** 

c 

c****** GENERATING NODES 1-121 
C 

11=0 

DO 60 1=1,11 
IA=1+ ( 1-1) *11 
YA=YY ( I ) 

XA=XX ( I ) 

AZ=XA* .15 
BZ=XA* .30 
CZ=XA* . 50 
DZ=XA* . 75 
XN ( IA) =-XA 
XN (IA+10) =XA 
XN ( IA+1 ) =-XA+AZ 
XN ( IA+9 ) =XA-AZ 
XN ( IA+2 ) =-XA+BZ 
XN ( IA+8 ) =XA-BZ 
XN ( IA+3 ) =-XA+CZ 
XN ( IA+7 ) =XA-CZ 
XN ( IA+4 ) =-XA+DZ 
XN ( IA+6 ) =XA-DZ 
XN (IA+5) =0 . 0 
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DO 60 J=l, 11 
11=11+1 
YN(II) =YA 
60 CONTINUE 
C 

C****** GENERATING NODES 122-209 
C 

DANG=2 . 5 
SANG=10. 0 
DO 70 1=1,8 
DXANG=SANG 
DO 75 J=l, 4 
XANG ( I , J ) =DXANG 
DXANG= DX ANG - D ANG 
WRITE (6 , *) XANG (I , J) 

XANG (I, J)=XANG(I,J)*PI/180. 

75 CONTINUE 

DANG=DANG+2 . 5 
SANG=SANG+10 . 

70 CONTINUE 

DO 71 1=1,8 
XANG (I, 4) =0.0 

71 CONTINUE 
C 

DO 80 1=1,8 

IA=122+(I-1) *11 
XA=XX ( 11+1 ) 

YA=YY ( 11+1 ) 

AZ=XA* . 15 
BZ=XA* .30 
CZ=XA* . 50 
DZ=XA* .75 
DZ=DZ* (1. + (I*.02) ) 

XN(IA)=-XA 
XN (IA+10) =XA 

XN ( IA+ 1 ) =-XA+AZ *COS ( XANG (1,1)) 

XN ( IA+ 9 ) =-XN ( IA+ 1 ) 

XN ( I A+ 2 ) =XN ( I A+ 1 ) + ( BZ - AZ ) *COS (XANG ( 1 , 2 ) ) 
XN ( IA+8 ) =-XN ( IA+2 ) 

XN ( IA+ 3 ) =XN ( IA+2 ) + ( CZ-BZ ) *COS (XANG (1,3)) 
XN ( IA+7 ) =-XN ( IA+3 ) 

XN ( IA+4 ) =XN ( IA+3 ) + ( DZ-CZ ) *COS (XANG (I , 4 ) ) 
XN ( I A+ 6 ) =-XN ( I A+4 ) 

XN ( IA+5) =0 . 0 
IF ( I . EQ . 7 ) CZ=CZ*1.1 
IF ( I . EQ . 8 ) CZ=CZ*1 . 15 
YN(IA) =YA 
YN (IA+10) =YA 

YN ( I A+ 1 ) = YA- AZ * S IN ( XANG (1,1) ) 

YN ( IA+9 ) =YN ( IA+1 ) 

YN ( I A+ 2 ) = YN ( I A+ 1 ) - ( BZ - AZ ) *SIN (XANG ( 1 , 2 ) ) 
YN ( IA+8 )=YN( IA+2) 

YN(IA+3)=YN(IA+2) -(CZ-BZ) *SIN (XANG ( I , 3 ) ) 
YN ( IA+7 )=YN( IA+3) 

YN ( IA+4 ) =YN ( IA+3 ) - ( DZ-CZ ) * SIN (XANG (1,4)) 
YN (IA+6) =YN (IA+4) 

YN ( IA+5 )=YN( IA+4) 

80 CONTINUE 
C 
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DZ.=DZ/1 = 16 
CZ=CZ/1 . 15 
C 

C****** GENERATING NODES 210-233 
C 

DO 85 1=1,3 
J=20+(I-1) 

XINC (I) =ATAN (XX ( J) /YY ( J) ) 

85 CONTINUE 

DINC=XINC(l)/2. 

DO 90 1=1,3 

IA=210+ (1-1) *8 
YA=YY ( 19+1) 

XA=XX(19+I) 

XN ( IA) =-XA 
XN ( IA+7 ) =XA 

XN (IA+1) =XN (IA) +AZ*SIN (XINC(I) ) 

XN ( I A+ 6 ) =-XN (IA+1) 

XN (IA+2 ) =XN (IA+1) + (BZ-AZ) *SIN(XINC(I) ) 

XN ( IA+5 ) =-XN ( IA+2 ) 

XN (IA+3) =XN (IA+2) + (CZ-BZ) *SIN(XINC(I) ) 

XN ( IA+4 ) =-XN ( IA+3 ) 

YN (IA) =YA 
YN ( IA+7 ) =YA 

YN ( IA+1) =YA-AZ*COS (XINC (I) ) 

YN (IA+6) =YN (IA+1) 

YN (IA+2 )=YN (IA+1) -(BZ-AZ) *COSfXINCm ) 

YN ( IA+5 j=YN( IA+2) 

YN ( IA+3 )=YN( IA+2) -(CZ-BZ) *C0S(XINC(I) ) 

YN ( IA+4 ) =YN ( IA+3 ) 

90 CONTINUE 
C 

C****** GENERATING NODES 234-241 
C 

DELTA=ATAN ( XX (23) / YY (23) ) 

ALPHA=DELTA-ATAN (XN(206)/YN(206)) 

XN (241) =XX (23) 

YN (241) =Y Y (23) 

XN (234) =-XN (241) 

YN (234) = YN (241) 

XN (240) =XN (241) -AZ*SIN (DELTA) 

YN (240) =YN ( 241) -AZ* COS (DELTA) 

XN (235) =-XN (240) 

YN (235) = YN (240) 

XN (239) =XN (241) -BZ*SIN (DELTA) 

YN ( 2 3 9 ) =YN ( 2 4 1 ) -BZ *COS ( DELTA) 

XN (236) =-XN (239) 

YN (236) =YN (239) 

XN (238) =XN (241) -CZ*SIN (DELTA) 

YN (238) = YN (241) -CZ*COS (DELTA) 

XN (237) =-XN (238) 

YN (237) = YN (238) 

C 

C****** GENERATING NODES 242-319 \ 

C 

RIM= (YN (238) -YY (27) ) /COS (DELTA) 

DRIM=RIM/6 . 0 
DO 100 1=1,5 
IA=238- (I — 1 ) *8 
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XADA(I) =ATAN (XN (IA)/YN (IA) ) 

XADA(6) =ATAN (XN (205)/YN (205) ) 

100 CONTINUE 
C 

DO 110 1=1,6 
IA=238— (1-1) *8 
IX=254- ( 1-1) 

IF (I . EQ. 6) IA=205 

XN (254) =XN (238) -DRIM*SIN (XADA(l) ) 

YN ( 2 5 4 ) =YN (238) -DRIM*COS ( XADA ( 1 ) ) 
XRAD=XN (254) /SIN (XADA ( 1) ) 

IF(I.EQ.l) GO TO 111 
XN (IX) =XRAD*SIN (XADA (I) ) 

YN (IX) =XRAD*COS (XADA(I) ) 

111 XN ( IX- (14-2*1) ) =-XN (IX) 

YN ( IX- (14-2*1) ) =YN (IX) 

110 CONTINUE 
C 

DO 119 J=1 , 6 
DO 120 1=1,6 
IA=255+(J-1) *13+ ( 1-1) 

XN ( IA) =XN ( IA-13 ) +DRIM*SIN ( XADA (I) ) 

YN ( I A) = YN ( I A- 1 3 ) -DRIM* COS ( XADA (I) ) 

XN (IA+ ( 14-2*1 ) ) =-XN (IA) 

YN (IA+ (14-2*1) ) =YN (IA) 

120 CONTINUE 

119 CONTINUE 

DO 125 1=1,6 

IA=248+ ( ( I — 1 ) *13) 

XN (IA) =0 . 0 
YN (IA) =YN (IA-1) 

125 CONTINUE 
C 

DO 130 1=1,319 

WRITE (11, 1000) I , XN ( I ) , YN ( I ) 

130 CONTINUE 

1000 FORMAT (I 10, 7X, '0.0' , F10 . 5 , F10 . 5) 

C 

C****** THIS ENDS THE NODE GENERATION SECTION 

C 

C 

C ***************************************** 

C * THIS SECTION GENERATES THE ELEMENTS * 

C * AND ELEMENT NUMBERS. TWO-DIMENSIONAL* 

C * 4-NODED ELEMENTS ARE USED. FOR THIS * 

C * MODEL THERE ARE 276 ELEMENTS. * 

C ***************************************** 

c 

c****** GENERATING ELEMENTS 1-183 
C 

11=0 

DO 140 1=1,18 
IA=12+ (I — 1 ) *11 
IB=13+ (1-1) *11 
IC=2+ ( I — 1 ) *11 
ID=1+ (1-1) *11 
DO 145 J=l,10 
11=11+1 
LA=IA+ ( J-l ) 
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LC=IC+ (J-l) 

LD=ID+(J-1) 

145 WRITE ( 11 / 112 ) II, LA , LB , LC , LD 

140 CONTINUE 

WRITE (11, 112) 181,210,211,200,199 
WRITE (11, 112) 182,211,212,201,200 
WRITE (11, 112) 183,212,213,202,201 
C 

C****** GENERATING ELEMENTS 184-195 
C 

11=183 

DO 150 1=1,3 
IA=218+ (1-1) *8 
IB=219+ (1-1) *8 
IC=211+ (1-1) *8 
ID=210+ (1-1) *8 
DO 155 J=l, 3 
11 = 11+1 
LA=IA+ (J-l) 

LB=IB+ (J-l) 

LC=IC+ (J-l) 

LD=ID+ (J-l) 

155 WRITE (11, 112) II, LA, LB, LC,LD 

150 CONTINUE 

WRITE (11, 112) 193,214,215,207,206 
WRITE (11, 112) 194,215,216,208,207 
WRITE (11, 112) 195,216,217,209,208 
C 

C****** GENERATING ELEMENTS 196-203 
C 

11=195 

DO 160 1=1,3 
IA=222+ (1-1) *8 
IB=223+ ( 1-1) *8 
IC=215+ (1-1) *8 
ID=214+ (1-1) *8 
DO 165 J=1 , 3 
11 = 11+1 
LA=IA+ (J-l) 

LB=IB+ (J-l) 

LC=IC+ (J-l) 

LD=ID+ (J-l) 

165 WRITE (11, 112) II, LA, LB, LC,LD 

160 CONTINUE 
C 

C****** GENERATING ELEMENTS 204-213 
C 

11=204 

DO 170 1=1,3 
LA=242+ (1-1) 

LB=243+ (1-1) 

LC=229- (1-1) *8 
LD=237- (1-1) *8 
11 = 11+1 

170 WRITE (11, 112) II, LA, LB, LC,LD 

WRITE (11, 112) 208, 245, 246, 202, 213 
WRITE (11, 112) 209,246,247,203,202 
WRITE (11, 112) 210,247,248,204,203 
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nniiti (ii, XJ.2 ) 2 J. 1 , 2 48 , 24 y , 2 05 , 204 
WRITE (11, 112)2 12 ,249, 250,206,205 
WRITE (11, 112) 213 ,250 ,251, 2 14, 206 
C 

C****** GENERATING ELEMENTS 214-276 
C 

11=213 

DO 180 1=1,3 
LA=251+(I-1) 

LB=252+ (1-1) 

LC=222+ (1-1) *8 
LD=2 14+ ( 1-1) *8 
11 = 11+1 

180 WRITE(11, 112)11, LA, LB, LC,LD 
11=216 

DO 190 1=1,5 

IA=255+ ( 1-1) *13 
IB=256+ ( 1-1) *13 
IC=243+ (1-1) *13 
ID=242+ ( 1-1) *13 
DO 195 J=l, 12 
11 = 11+1 
LA=IA+ (J-l) 

LB=IB+ ( J-l) 

LC=IC+ (J-l) 

LD=ID+ (J-l) 

195 WRITE (11, 112) II, LA, LB, LC,LD 

190 CONTINUE 
C 

112 FORMAT (15, 4X, '8' ,4I5,23X, *21' ,4X, '1') 

C 

c****** THIS ENDS THE ELEMENT GENERATION SECTION 
C 

STOP 

DEBUG SUBCHK 
END 
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onoooooooooooooooooooooooononooooon 


C*** ************************************** ********* 


* 

VARIABLE SPEED/LOAD GENERATION PROGRAM * 

* 

THIS PROGRAM GENERATES THE SAPS TIME * 

FUNCTIONS DESCRIBING THE LOAD MOVING * 

OVER THE INVOLUTE OF A GEAR TOOTH. THE * 

SPEED OF THE MOVING LOAD ON THE INVOLUTE * 

OF A TOOTH VARIES AS * 

* 

V (T) =RB*OMEGA**2*T * 

* 

WHERE RB IS THE BASE RADIUS OF THE GEAR * 

OMEGA IS THE ANGULAR VELOCITY OF THE GEAR * 

AND T IS THE ELLAPSED TIME. THE LOAD * 

FUNCTIONS WERE TAKEN FROM WALLACE, AND ARE * 

* 

F (T) =F0 ( 1-COS (T/ALPHAT) ) /2 0<T<ALPHAT * 

* 

F (T) =F0 ALPHAT<T<( 1-ALPHA) TF * 

* 

F (T) =F0 ( 1-COS (T-ALPHATF) /ALPHATF) )/2 * 

( 1-ALPHA) TF<T<TF * 

* 

WHERE TF IS THE TOTAL TIME, AND ALPHA IS * 

A FACTOR DEPENDENT ON THE CONTACT RATIO. * 

* 

THE USER MUST ENTER THE FOLLOWING PARAMETERS * 

* 

PRESSURE ANGLE=PANG * 

PITCH RADIUS =RP * 

ADDENDEM =AD * 

CIRCULAR PITCH=CIRP * 

BACKLASH =BACKL * 

MAX VELOCITY =VMAX * 

* 


C* ************************************************* 

DIMENSION XX (40) ,YY(40) , FORCE (11 , 500) ,TIME(500) ,T(500) 

DIMENSION XINC (11) ,DXINC(11) ,FX(11,500) ,FY(11,500) ,XN(15) ,YN(15) 
PI=3. 141592654 
READ ( 5 , * ) PANG 
PANG=PANG*PI/180 . 

READ (5,*) RP 
READ ( 5 , * ) AD 
READ (5,*) CIRP 
READ (5,*) BACKL 
READ (5,*) VMAX 
RB=RP*COS (PANG) 

THP= ( (RP*RP-RB*RB) / (RB*RB) ) **0 . 5 
CIRTH=CIRP/2 . -BACKL 
AL=THP-ATAN (THP) +CIRTH/ (RP*2 . ) 

RO=RP+AD 

C 

C* ************************************************* 

C THE INVOLUTE COORDINATES AND INVOLUTE NORMALS * 

C ARE CALCULATED IN THIS SECTION. * 

C********** ****** ******* ******************** ******* 


DO 5 1=1,11 
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IA=1+ (1-1) 

HH= (RO-RB) /10 . 

H=RO-HH* (1-1) 

ANG= ( (H*H-RB*RB) / (RB*RB) ) **0.5 
Y=RB* (COS (ANG) +ANG*SIN (ANG) -1 . ) 

X=-RB* (SIN (ANG) -ANG*COS (ANG) ) 

XP=RB*SIN (AL) +X*COS ( AL) +Y*SIN (AL) 

YP=RB*COS (AL) -X*SIN ( AL) +Y*COS ( AL) 

XX (IA) =XP 
YY (IA) =YP 

XN (I) =-SIN (AL) *SIN (ANG) -COS (AL) *COS (ANG) 

YN (I) =SIN (AL) *COS (ANG) -COS (AL) *SIN(ANG) 

5 CONTINUE 
C 

C************************************************** 

C THE DISTANCE ALONG THE INVOLUTE IS DETERMINED * 

C FROM WHICH THE TOTAL TIME AND ANGULAR VELOCITY * 

C ARE FOUND. * 

C************************************************** 

C 

S=0. 0 

DO 6 1=1,10 

XL=XX ( 1+1) -XX (I) 

YL=YY (I) -YY (1+1) 

XINC ( I ) = (XL**2+YL**2 ) **0.5 
S=S+XINC (I) 

WRITE ( 6 , * ) ' XINC (I) = ' ,XINC(I),' S= ',S 

6 CONTINUE 
TF=2 . *S/VMAX 

XOMEGA= (VMAX/ (TF*RB) ) ** . 5 

WRITE (6,*) 1 TF= ' ,TF, 1 OMEGA= ' , XOMEGA 

TALPHA=0 . 2 8 *TF 

WRITE (6,*) ' TALPHA= ' , TALPHA , * ( 1-ALPHA) *TF= ' , ( 1 . - . 28) *TF 

C 

C* ************************************************* 

C THE DISTANCE BETWEEN EACH NODE IS DIVIDED UP * 

C INTO 20 EQUAL INCREMENTS AND THE TIME AND * 

C AND FORCE VALUES ARE DETERMINED FOR EACH. * 

C**************** ******************* *************** 

C 

DIST=0 . 0 
DO 10 1=1,10 

DXINC (I) =XINC (I ) /20 . 

DO 15 J=1 ,20 
JJ=J+20*(I-1) 

DIST=DIST+ DXINC ( I ) 

TIME ( JJ) = (2 . *DIST/ (RB* XOMEGA* *2) ) **0 . 5 
T(JJ)=TIME(JJ) 

IF (TIME (JJ) .GT. TALPHA) GO TO 11 

FORCE (I , J) = ( FO/2 . ) * ( 1 . -COS ( PI *TIME ( J J ) /TALPHA) ) 

GO TO 15 

11 IF(TIME(JJ) .GT. (TF-TALPHA) ) GO TO 12 
FORCE (I, J)=FO 

GO TO 15 

12 FORCE ( I, J)= (FO/2. ) * (1. -COS (PI* (TF-TIME (JJ) ) /TALPHA) ) 

15 CONTINUE 

10 CONTINUE 
C 


DO 20 1=1,1 
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DO 30 J=l, 20 
JJ=J+20* (1-1) 

FX (I , J+20) =FORCE (I , J) * ( (XINC (I) - (RB* (X0MEGA**2 ) /2 . ) 

$*TIME(JJ) **2) /XINC(I) ) *XN (I) 

FY (I , J+20) =FORCE (I , J) * ( (XINC (I) - (RB* (XOMEGA**2 ) /2 . ) 

$*TIME(JJ) **2 ) /XINC (I) ) *YN (I) 

FACT=(RB*(XOMEGA**2)/2. ) * (TIME ( JJ+180) **2-TIME ( 180) **2) 
FX(11, J)=FORCE(10, J) * (FACT/ XINC (10) ) *XN(11) 

FY ( 11 , J) =FORCE ( 10 ; J) * (FACT/XINC (10) ) *YN(11) 

30 CONTINUE 

20 CONTINUE 

DO 25 1=2,10 
DO 35 J=1 ,20 
JJ=J+20* (1-2 ) 

XTIME1=TIME ( (1-1) *20) 

IF (I . NE . 2 ) GO TO 36 
XTIME=0 . 0 
GO TO 37 

36 XTIME=TIME( (1-2) *20) 

37 FACT=( (RB*(XOMEGA**2) )/2. )* (TIME(JJ) **2-XTIME**2 ) 

FACT1=( (RB*(XOMEGA**2) )/2.) * (TIME ( JJ+20) **2-XTIMEl**2 ) 

FX (I , J) =FORCE (1-1 , J) * (FACT/XINC (1-1) ) *XN (I) 

FY (I, J)=FORCE(I-l,J) * (FACT/XINC (1-1) ) *YN (I) 

FX ( I , J+2 0 ) =FORCE (I , J) * ( ( XINC ( I ) -FACT1) /XINC (I) ) *XN (I) 

FY ( I , J+2 0 ) =FORCE ( I , J) * ( ( XINC ( I ) -FACT1) /XINC ( I ) ) * YN ( I ) 

35 CONTINUE 

2 5 CONTINUE 
C 

C************************************************* 

C THE REMAINDER OF THE PROGRAM WRITES THE TIME * 

C FUNCTION VALUES (TIME, FORCE) TO A FILE FOR * 

C USE IN A FINITE ELEMENT CODE. THIS PROGRAM * 

C IS FORMATTED FOR SAP6 . * 

C* ************************************ ************ 

C 

NFN=22 

WRITE (11, 99) NFN 

99 FORMAT (3X, 12, 4X, ' 0 ' , 4X, • 0 ' , 3X, ' NT • , 4X, • 1 ' , 6X, 'DPER',7X, '0.0') 

IFN=0 

DO 60 1=1,11 
NP=11+11* (1-1) 

IC=2 

DO 60 J=1 , 2 
IFN=IFN+1 

WRITE (11,98) NP, IC, IFN 
IC=3 

98 FORMAT (15, 4X, II, 15, 4X, '1\7X, '1.0',24X, '0') 

60 CONTINUE 

WRITE (11, 97) 

97 FORMAT (/' ') 

C 

C***** TIME FUNTION DATA **** 

C 

TF=1 . 0 
T0=0 . 0 
F0=0 . 0 

C******** TIME FUNCTION FOR NODE 11 ******* 

C 


NLP=24 
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WRITE (11, 101) NLP 

101 FORMAT (3X, 12, 11X, 'TIME FUNCTION ', 50X, ' 1 ') 

100 FORMAT (3 (E12 . 6, F12 • 0) ) 

WRITE ( 11 , 100) TO , FO , T ( 1) ,FX(1,21) ,T(2) ,FX(1,22) 

DO 200 1=3,18,3 

200 WRITE (11,100) T (I) ,FX(l,I+20) ,T(I+1) ,FX(1,I+21) ,T(I+2) ,FX(l,I+22) 
WRITE (11,100) TF , FO , TF , FO , TF , FO 

WRITE (11, 101) NLP 

WRITE (11,100) TO , FO , T ( 1) ,FY(1,21) ,T(2) ,FY(1,22) 

DO 201 1=3,18,3 

201 WRITE ( 11 , 100) T (I) ,FY(l,I+20) ,T(I+1) ,FY(1,I+21) ,T(I+2) ,FY(l,I+22) 
WRITE ( 11 , 100) TF , FO , TF , FO , TF , FO 

C****** TIME FUNCTIONS FOR NODES 22 THRU 110 ***** 

C 

DO 202 1=2,10 
NLP=45 

WRITE (11, 101) NLP 
K=20* (1-2 ) 

T ( 0) =0 . 0 

WRITE (11, 100) TO , FO , T (K) ,F0,T(1+K) ,FX(I,1) 

DO 203 J=2 ,38,3 

203 WRITE ( 11 , 100) T(J+K) ,FX(I,J) ,T(J+1+K) ,FX(I,J+1) ,T(J+2+K) ,FX(I,J+2; 
WRITE (11,100) TF, FO , TF, FO ,TF, FO 
WRITE (11, 101) NLP 

WRITE (11, 100) TO , FO , T (K) ,F0,T(1+K) ,FY(I, 1) 

DO 205 J=2 ,38,3 

205 WRITE ( 11 , 100) T(J+K) ,FY(I,J) ,T(J+1+K) ,FY(I,J+1) ,T(J+2+K) ,FY(I,J+2; 
WRITE ( 11 , 100) TF, FO , TF, FO, TF, FO 

202 CONTINUE 
C 

c**** TIME FUNCTION FOR NODE 121 ******* 

c 

NLP=24 

WRITE (11, 101) NLP 

WRITE (11,100) TO , FO , T ( 180) ,F0,T(181) ,FX(11,1) 

DO 206 1=2,17,3 
K=I+180 

206 WRITE (11, 100) T (K) , FX ( 11 , I) ,T(K+1) ,FX( 11,1+1) ,T(K+2) ,FX(ll,I+2) 
WRITE (11,100) T (200) ,FX(11,20) ,TF,F0,TF,F0 

WRITE (11, 101) NLP 

WRITE (11,100) TO, F0,T(180) ,F0,T(181) ,FY(11,1) 

DO 207 1=2,17,3 
K=I+180 

207 WRITE (11, 100) T (K) , FY ( 11 , 1) ,T(K+1) , FY (11,1+1) ,T(K+2) ,FY(ll,I+2) 
WRITE (11, 100) T (200) ,FY(11,20) ,TF,F0,TF,F0 

C 

C**** ECHO CHECK ****** 

C 

DO 40 1=1,11 

IF(I.NE.l) GO TO 52 
DO 42 J=1 , 20 
JJ=J+20*(I-1) 

WRITE (12,51) T(JJ) ,FX(I,J+20) 

42 CONTINUE 

GO TO 40 

52 IF(I.EQ.ll) GO TO 53 

DO 45 J=1 , 20 
JJ=J+20*(I-2) 

WRITE ( 12 , 51) T(JJ) ,FX(I, J) 
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45 

46 
53 

50 
40 

51 


CONTINUE 
DO 46 J=l, 20 
JJ=J+20*(I-1) 

WRITE (12,51) T(JJ) ,FX(I, J+20) 
CONTINUE 
GO TO 40 
DO 50 J=1 , 2 0 
JJ=J+20* (1-2) 

WRITE (12,51) T(JJ) ,FX(I, J) 
CONTINUE 
CONTINUE 

FORMAT (F12 . 8 , 12X, F8 . 0) 

STOP 

END 
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APPENDIX 2 

Contact Point Velocity 
Of Meshing Spur Gear Teeth 
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For a meshing spur gear pair, the speed of the contact 
point varies with time as it moves from the tip of the tooth to 
the base circle. 

From equations (2-3) and (2-4) and Figure A2-1 the coor- 
dinates of a point Bi on the involute are given by; 

XBi * -RB( sin 0i - 9iCOS9i) (A2-1) 

YBi = RB( COS0 x + 0isin0i - 1) (A2-2 ) 


Differentiating (A2-1) and (A2-2) with respect to 0j.; 

dXB . „ „ . A 

i= - RB0iSin0i 


dYB ± 

a ~ Q ~ - RB0icos0i 


and combining to give the resultant; 


9 dXB, 9 dYB. 9 
/ dr '2_, l , 2 , i.2 

'd0 . 1 d0 . 1 d0 . ' 

ii i 


= RB 2 ( sin 2 0i+cos 2 0i) 


(A2-3 ) 
(A2-4 ) 


= RB0i (A2-5 ) 

Multiplying through by d 0i = 0j_, realizing that 0f and r vary 
with time, we have 

= V = RB0i9i 

And if we define 0 i = cot; 



where; 


V(t) = RBoj 2 t 


(A2-6 ) 
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Figure A2-1: Contact point goemetry 
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V(t) = the speed of the contact 
point Bi at time t 


w = angular velocity of the gear (rad/sec) 

At t = 0.0, corresponding to the load at the tip of the tooth, 
the speed is zero, and at t = TF, time corresponding to the 
load positioned at the base circle, the speed is maximum. 

The position of the load as a function of time is found by 
integrating equation (A2-6) with respect to time. By defining 
the speed as the first derivative of the position; 

V( t) = dSjt) _ RBu) 2 t 
dt 

o f 

dS ( t ) = RB of j t dt 

S(t) = ^RBoj2t2 + C ( A2-7 ) 


Evaluating the constant of integration for t ■ Tj., the time at 
the tip of the tooth; 


then; 


S(ti) = 0.0 = ^RBco2T i 2+ c 


C = -^RBw 2 t| 

The displacement then becomes; 

S(t) = J£RBa) 2 t2 - ^2RBo)2Ti2 


(A2-8 ) 


(A2-9 ) 


If we assume that the time at the tip position is zero the 
constant term vanishes and we are left with only the first term 
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on the right hand side of equation A2-9 defining the position 
along the involute. 
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APPENDIX 3 

Deformed Shapes of Gear Tooth 

1. Static Loading 

2. Modal Analysis 

3 . Dynamic Response 
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Static Loading 
Timoshenko Beam Constraints 
Rim Included 
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STRTIC DEFLECTION ANALYSIS : GEAR »1 (BEAM CONSTR) 
STATIC LOAD CASE 1 t = O.l 


SEPTEMBER 09. 1983 01:37:07 

IRXIS=3 RLPHfl= 0.00 BETR= 0.00 

DEFLECTION SCALE FRCTQRKI + . 1 


z 


% 


Y 



Figure A3-2 
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STATIC DEFLECTION ANALYSIS: GEAR #1 (BEAM CONSTR) 
STATIC LOAD CASE 2 t = 0.3 


SEPTEMBER 09. 1983 01:37:07 

IRX1 S=3 RLPHR= 0.00 BETfl= 
DEFLECT I ON SCALE FACTOR 8 .6 


■ 


M 


■n 


■■■ 


Figure A3-3 


STATIC DEFLECTION ANALYSIS: GEAR #1 (BEAM CONSTR) 
STATIC LOAD CASE 3 t = 0.5 

SEPTEMBER 09. 1983 01:37:07 

IAXIS=3 ALPHfi= 0.00 BETR= 0.00 
DEFLECTION SCALE FRCTCK&67 .8 



Figure A3 -4 
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STATIC DEFLECTION ANALYSIS: GEAR #1 
STATIC LOAD CASE 4 T = 0.8 

SEPTEMBER 09 . 1983 0 1 *37 :07 

IHXIS=3 ALPHA= 0.00 BETA= 0.00 

DEFLECTION SCALE FACTO®? . 


(BEAM CONSTR ) 




HobB 

■^n 


Figure A-5 
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STATIC DEFLECTION ANALYSIS: GEAR #1 
STATIC LOAD CASE 1 t = o.l 

SEPTEMBER 10, 1983 00:10:05 

IFIXIS=3 BLPHB= 0.00 BETfl= 0.00 

DEFLECTION SCALE FACTORS 3 .8 


(RIM INCLUDED) 

2 

X Y 



Figure A3 -6 
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STATIC DEFLECTION ANALYSIS: GEAR #1 (RIM INCLUDED) 
STATIC LOAD CASE 4 t = 0.8 


EPTEMBER 10. 1983 00:10:05 

HX1 S=3 RLPHFI= 0.00 BETH= 0.00 

EFLECTION SCALE FRCTCKWH .3 


% 


Y 



Figure A3-9 
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Modal Analysis 
Timoshenko Beam Constraints 
Rim Constraints 
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MODAL ANALYSIS, PLANE STRAIN ELEMENT: GEAR #1 
UNDEFORMED SHAPE 


UNE 1G. 1983 00:10:15" 

1RXIS=3 HLPHfl= 0.00 BETfl= 0.00 



Figure A3-10' 
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MODRL ANALYSIS, PLANE STRAIN ELEMENT: GEAR #1 
MODE 1 Freq.= .9101 E5 (CPS) 

UNE 10, 1383 13:32:39' 

inXIS=3 RLPHR= 0.00 BETR= 0.00 

DEFLECTION SCALE FRCTOR= 0.24278-E 



Figure A3-11 
2 


MODAL ANALYSIS, PLANE STRAIN ELEMENT: GEAR #1 
MODE 2 Freq.= .2159 E6 (CPS) 

UNE 10. 1983 19:32:39' 

!RXIS=3 RLPHA= 0 .OH BETH- 0.00 
DEFLECTION SCALE FRCTOR= 0 . 30998-E 



Figure A3-12 
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MODRL RNRLYSIS, PLRNE STRRIN ELEMENT: GERR # 
MODE 3 Freq.= .2289 E6 (CPS) 

UNE 10, 1983 1 9 5 3 2 : 3 9 * 

1 PXI S'=3 nLPHH= o.on BEI'R= 0.00 
DEFLECTION SCALE FACTOR= 0 . 24 7 5 6-E 



MODAL ANALYSIS. PLANE STRAIN ELEMENT: GEAR #1 
MODE 4- Freq.= .4218 E6 (CPS) 

UNE 10. 1983 19:32:39' 

1 HXI S=3 RL.PHR= 0.00 BET R= 0.00 

DEFLECTION SCHLE FHCTOR= 0 . 2 1 H 4 4-E 



Figure A3-14 
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MODAL ANALYSIS, PLANE STRAIN ELEMENT: GEAR #1 
MODE 5 Freq.= .4794 E6 (CPS) 



Figure A3-15 
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MODRL ANALYSIS, PLANE STRAIN ELEMENT 
MODE 1 Freq.= .6545 E5 (CPS) 


GEAR #1 


UNE 16, 1983 00:10:15' 

I AXI S=3 ALPHA= 0.0 0 BETA= 
DEFLECT I ON SCALE FACTOR^ 0 . 27565-E 


OH 


Figure A3-16 



MODAL ANALYSIS. PLANE STRAIN ELEMENT: GEAR #1 
MODE 2 Freq . = .1146 E6, (CPS) 

UNE 16. 1983 00 s 10 : 15 1 

IRXIS=3 RLPHfl= 0.00 BETB= 0.00 

DEFLECTION SCRLE FFICTOR= 0 .54125C 



Figure A3-17 
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MODAL ANALYSIS. PLANE STRAIN ELEMENT: GEAR #1 
MODE 3 Freq .= .1517. E6 (CPS) 


UNE IS, 1383 00:10:15' 

I AXI S=3 ALPHA= 0.00 BETA= 0.00 
DEFLECTION SCALE FACTOR= 0 .3237 EE 


z 


-¥ 



Figure A3-18 
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MODRL ANALYSIS, PLANE STRAIN ELEMENT 
MODE 1 Freq.= .2034 E6 (CPS) 

UNE 16. 1983 00 : 10 : 15 ' 

IRXIS=3 RLPHR= 0.00 BETR= 0.00 
DEFLECTION SCRLE FRCTOR= 0 . 1 2 67fE 


GEAR #1 


m 




Figure A3-19 


MODAL ANALYSIS, PLANE STRAIN ELEMENT: GEAR #1 
MODE 5 Freq.= .2319 E6 (CPS) 

UNE 16 . 1983 00 * 10 s 15 1 

IRXIS=3 RLPHR= 0.00 BETR= 0.00 

DEFLECTION SCALE FRCTOR= 0 . 5 2 5 3 8-E 



Figure A3-20 
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Dynamic Response 


Gear Tooth Profile Deflections 
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In the modal response analysis , modal superposition is 
used to sum the effects of all the included mode shapes pro- 
ducing the desired response. The more modes included, the more 
accurate the analysis. However, there is a practical limit to 
the number of modes used dictated by the type of loading used 
(impact, constant, sinusoidal, etc.) and of course com- 
putational efficiency. As previously mentioned, the first ten 
modes are used for this analysis. 

To better comprehend the dynamic deflection phenomena, the 
deflected state of the tooth profile, subject to the impact 
loading case, is plotted for different load positions with 
V*=0.01. Figure A3-21 shows the underformed shape and Figures 
A3-22 through A3-25 are the various deformed shapes. From the 
figure(s), one can see those modes which contribute noticeably 
to the overall deflection of the tooth. It is apparent that 
only the first three or four have a major effect (see Figure 
A3-19 , T=0.5). Thus the problem of local deformation encoun- 
tered in the static analysis is not apparent here. 
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DYNAMIC DEFLECTIONS OF GEAR TOOTH: UNDEFORMED SHAPE 


Figure A3-21 




V*=0.01 T=0.1 

DYNAMIC DEFLECTIONS OF GEAR TOOTH WITH MOVING LOADS 


Figure A3-22 
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v*=0.01 


T=0.3 


DYNAMIC DEFLECTIONS OF GEAR TOOTH WITH MOVING LOADS 


Figure A3- 2 3 
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V*=0.01 T=0.5 

DYNAMIC DEFLECTIONS OF GEAR TOOTH WITH MOVING LOADS 


Figure A3-24 
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V*=0.01 


T=o.a 


DYNAMIC DEFLECTIONS OF GEAR TOOTH WITH MOVING LOADS 


Figure A3-25 
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DISPLACEMENT(INCHES) DISPLACEMENT(INCHES) DISPLACEMENT (I NCHES) 
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POSITION 

VEL= 10.0 IN/SEC ; MASS= 1.0 ; BEAM LENGTH^ 0.5 IN 
DYNAMIC RESPONSE OF MESHING CANTILEVER BEAMS, WITH INERTIA (VEL=V(T)) 


Figure A4-16 
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DISPLACEMENT(INCHES) DISPLACEMENT NCHES) DISPLACEMENT(INCHES) 
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VEL= 20.0 IN/SEC ; MASS= 1.0 ; BEAM LENGTH = 0.5 IN 


DYNAMIC RESPONSE OF MESHING CANTILEVER BEAMS, WITH INERTIA (VEL=V(T)) 


Figure A4-17 
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DISPLACEMENT(INCHES) DISPLACEMENT(INCHES) DISPLACEMENT(INCHES) 


.25 


.20 ■ 
.15 - 





VEL= 40.0 IN/SEC ; MASS= 1.0 ; BEAM LENGTH— 0.5 IN 


DYNAMIC RESPONSE OF MESHING CANTILEVER BEAMS. WITH INERTIA (VEL=V(T)) 


Figure A4-18 











APPENDIX 5 


Use of Natural Modes in Equation 
of Motion for Cantilever Beams 

1) Constant of Integration 

2) Evaluation of I2PHI1 AND I2PHI2 

3) Derivatives of the Mode Shape With Time 



1.) Constant of Integration 

In order to simplify the solution of the equation of 
motion (5-26), the constant of integration associated with the 
mode shape of a cantilever beam is determined such that; 

p( L <)> 2 (Ud£ = 1 (A5-1) 

)o 

Since the beams are identical and the integral is over the 
length of the beam, the constant is the same for both beams. 

The first natural mode of vibration for a uniform fixed- 
free cantilever beam is; 


<{>(?) = C[ (sin$L-sinh3L) ( sin $£-sinh 35) 


+ ( cos gL+cosh gL ) ( cos 35 -cosh g£) 


(A5-2 ) 


where; 


4 - 

e = 


2 

CO p 

El 


Equation (A5-2) is simplified by rewriting with; 

51 = singL-sinhgL 

52 = cosgL+coshfSL 


We then have; 

<J)( £) = C[ Slsinx-Slsinhx+S2cosx-S2coshx] (A5-3) 

where x is used in place of the argument g£ . 
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Squaring equation (A5-3); 


<t> 2 U) = c2 [Sl2 ( sin2x+sinh2x-2sinxsinhx) 

+ S22 ( cos2x+cosh2x-2cosxsinhx) (A5-4) 

+2SlS2(sinx cosx-sinx coshx-cosx sinhx+sinhx coshx)] 

We can then solve for the constant C with the following 
substitution; 


C =(— — ) h (A5-5) 

P f G(C)d £ 

0 

where <j> ( E, )=CG( 5). Whenever the relation represented by 
equation (A5-1) occurs in equations (5-26) it is replaced by 1. 
Throughout the remainder of equations of motion the constant of 
integration C is represented by equation (A5-5). This relation 
is evaluted using Simpson's Rule. 

2.) Evaluation of I2PHI1 and I2PHI2 

Using the arbitrary constant of integration, C, determined 

in the previous section, the relationship; 

L 

I2PHI1 = 1 2 PH 1 2 =f <)> 2 (£)d£ (A5-6) 

J 0 

can be evaluated. The integrand represents the second deriva- 
tive of the mode shape with respect to the length variable, £ , 
quantity squared. The first derivative of the mode shape with 
x again used in place of 

= Cfl[Slcosx-S2coshx-S2sinx-S2sinhx] 

d£ r 
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and the second derivative; 


Squaring 

result; 



this 


= C£ 2 [-Slsinx-Slsinhx-S2cosx-S2coshx] 
relation and simplifying yields the 


desired 


2 9 

( d <jK 2 _ c 2g4 [Sl(sin 2 x+sinh z x+2sinx sinhx) 
+S2 2 ( cos 2 x+cosh 2 x+2cosx coshx) 
+2SlS2(sinx cosx+sinx coshx+cosx sinhx+sinhx coshx 


(A5-7) 

)] 


This relation is then used in equation (A5-6) which is eva- 
luated using Simpson's Rule. 


3.) Derivatives of the Mode Shape with Time 


When the constraint equation is differentiated, first and 


second derivatives of the mode shape result where the dependent 
variable, £ , must be considered as a function of time. The 
first derivative of the mode shape with respect to time is; 


djm 

dt 


C££ [Slcosx-Slcoshx-S2sinx-S2 sinhx] 


( A5-8 ) 


where x is used in place of the argument £5 and i is the velo- 
city of the moving load. Taking the second derivative, we 
have; 


d 2 (J>(0 


dt 


C[Sl(-£ 2 £ 2 (sinx+sinhx) -f^cosx-coshx) ) 

- S2 ( B 2 £ 2 ( cosx+cosh ) 

+ ££ ( sinx+sinhx ) ) ] (A5-9) 


where £ is the acceleration of the moving load. For a moving 
load with constant speed, £ is of course zero. 
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APPENDIX 6 


Development and Solution of 
Equations of Motion with 
Inertial Terms Removed 



By eliminating the inertial terms from the first four 
equations of (5-26), and using the undifferentiated form of the 
constraint equation (5-18), a new set of equations for the 
beam-mass system take the form; 


[A]{X} = {B} 


(A6-1) 


where; 


[A] = 


(M1+MB1) 

0 

0 

0 

0 


{X} = 


(M2+MB2 ) 
0 
0 
0 



0 

0 

I 2 PH II 
0 

-PHI1 


{B} = 


0 

0 

0 

I2PHI2 
-PHI 2 



1 

-1 

PHI1 
PHI 2 
0 


The initial conditions for this system are determined by 
arbitrarily chosing three of the four unknowns composing the 
constraint equation. Letting the beams assume static deflec- 
tions at time equal to zero, and XI equal to zero; the initial 
conditions are; 


ql ( 0 ) 


P«PHI1 

I2PHI2 


q2 ( 0 ) 


P*PHI2 
= I2PHI2 


with; 


X1(0) = 0.0 
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X2 ( 0 ) = XI + qlPHIl + q2PHI2 


Solving for {X} in equation (A6-1) by inversion; 

{X} = [Al-iiB} 

gives the values of the unknowns in terms of XI and X2 . Using 
these values for {X}={X1 X2 ql q2 A } T , the system is integrated 
using a Runge-Kutta integration algorithm solving for XI and 
X2 and X. Then equations 3 and 4 of (A6-1) are used to solve for ql 
and q2 . 

Tests were done for constant velocity moving loads of 1.0, 
5.0, 10.0, 20.0 and 40.0 inches per second. Plots of X2-X1, 
qlPHIl and q2PHI2 are included in Figures (A6-1) through 
( A6-5 ) . Comparing these to Figures (5-7) to (5-11) show that 
even though only one vibrational mode is used, the results are 
very nearly the same as those determined for the massless 
beams . 
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DISPLACEMENT^ NCHES) DISPLACEMENT(INCHES) DISPLACEMENT(INCHES) 



TIME(SEC) 



TIME(SEC) 



TIME(SEC) 

VEL= 1.0 IN/SEC ; MASS= 1.0 ; BEAM LENGTH = 0.5 IN 


DYNAMIC RESPONSE OF MESHING CANTILEVER BEAMS, NO INERTIA (VEL=CONST) 

Figure A6-1: 
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DISPLACEMENT(INCHES) DISPLACEMENT (I NCHES) DISPLACEMENT (I NCHES) 


.25 


.20 h 



-.05 1 ' 1 1 1 1 1 1 *— 

.000 .010 .020 .030 .040 .050 .060 .070 .080 

TIME(SEC) 




TIME(SEC) 

VEL= 5.0 IN/SEC ; MASS= 1.0 ; BEAM LENGTH= 0.5 IN 


DYNAMIC RESPONSE OF MESHING CANTILEVER BEAMS, NO INERTIA (VEL=CONST) 


Figure A6-2: 
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DISPLACEMENT(INCHES) DISPLACEMENT(INCHES) DISPLACEMENT(INCHES) 


.25 


.20 [ 



TIME(SEC) 



TIME(SEC) 



TIME(SEC) 

VEL= 10.0 IN/SEC ; MASS= 1.0 ; BEAM LENGTH= 0.5 IN 


DYNAMIC RESPONSE OF MESHING CANTILEVER BEAMS, NO INERTIA (VEL=CONST) 


Figure A6-3 
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DISPLACEMENT(INCHES) DISPLACEMENT(INCHES) DISPLACEMENT(tNCHES) 





TIME(SEC) 

VEL= 20.0 IN/SEC ; MASS= 1.0 ; BEAM LENGTH = 0.5 IN 


DYNAMIC RESPONSE OE MESHING CANTILEVER BEAMS, NO INERTIA (VEL=CONST) 

Figure A6-4: 
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DISPLACEMENT^ NCHES) DiSPLACEMENT(INCHES) DISPLACEMENT(INCHES) 





TIME(SEC) 

VEL= 40.0 IN/SEC ; MASS= 1.0 ; BEAM LENGTH= 0.5 IN 


DYNAMIC RESPONSE OF MESHING CANTILEVER BEAMS, NO INERTIA (VEL=CONST) 


Figure A6-5: 
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APPENDIX 7 

Dynamic Response Algorithm 
for Meshing Cantilever Beam 

1. Massless Configuration 

2. Inertia of Beam Included 
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c*********************************************** 

C THIS PROGRAM SOLVES THE EQUATION OF MOTION * 

C DESCRIBING A PAIR OF MESHING CANTILEVER * 

C BEAMS ATTACHED TO MOVEABLE FOUNDATAION * 

C MASSES. THE SECOND ORDER DIFFERENTIAL * 

C EQUATION: M ( D2X/DT2 ) +KX = P IS WRITTEN * 

C AS A FIRST ORDER DIFFEQ AND INTEGRATED * 

C USING A 4TH ORDER "RUNGE-KUTTA" INTEGRA- * 

C TION ROUTINE. * 

C* ********************************************** 

EXTERNAL FCT,OUTP 

DIMENSION Y (4 ) , DERY ( 4 ) , PRMT (5) ,AUX(8,4) 

REAL P , L , BA, VEL, Ml , M2 , MASS , INERT , K1 , K2 , STIFF , MOD , STIFF1 
COMMON P , L , VEL , MASS , INERT , K1 , K2 , STIFF , MOD , ZETA1 , ZETA2 , TF , DENOM 
WRITE (6,*) 'ENTER: APPLIED LOAD, BEAM LENGTH, BASE WIDTH' 

READ ( 5 , * ) P, L, BA 

WRITE (6,*) 'ENTER: VELOCITY, MASS#1, MASS #2' 

READ ( 5 , * ) VEL, Ml, M2 
MOD=30.E6 

INERT* (1./12,. ) *BA*BA**3 
ZETA1=L/10. 

ZETA2=9. *L/10. 

Kl= ( 3 . *MOD*INERT) / ( ZETA1**3 ) 

K2=(3.*MOD*INERT)/(ZETA2**3) 

STIFF* (K1*K2 ) / (K1+K2 ) 

STIFF1=STIFF 
WRITE ( 6 , * ) 

WRITE (6, *) 

WRITE (6,*) 'BASE 
'MASS 
'MASS 


■ ,u uut\L 

'BEAM LENGTH 
WIDTH 
#1 
#2 


/ * 

,L 
, BA 
/Ml 
,M2 
, VEL 


' ,PRMT(1) 
' , PRMT ( 2 ) 
• , PRMT ( 3 ) 


WRITE (6 , *) 

WRITE (6,*) 

WRITE (6,*) 'VELOCITY 
Y ( 1) =0 . 0 
Y(2)=P/STIFF 
DERY (1) =0.5 
DERY ( 2 ) =0 . 5 
PRMT (1) =0.0 
PRMT (2)=(.8*L) /VEL 
PRMT ( 3 ) =PRMT ( 2 ) /5000 . 

WRITE (6,*) 'STARTING TIME: 

WRITE (6,*) 'ENDING TIME: 

WRITE (6,*) 'TIME INCREMENT: 

PRMT(4) =0.0001 
NDIM=2 

MASS* (M1*M2 ) / (M1+M2 ) 

TIME=0 . 0 

CALL RKGS ( PRMT , Y , DERY , NDIM , IHLF , FCT , OUTP , AUX ) 
******************************************************** 
******************************************************** 

SUBROUTINE FCT (T, Y , DERY) 

DIMENSION Y (4), DERY (4) 

REAL P , L , BA , VEL, Ml , M2 , MASS , INERT , K1 , K2 , STI FF , MOD 

COMMON P , L, VEL, MASS , INERT , K1 , K2 , STIFF , MOD , ZETA1 , ZETA2 , TF , DENOM 

ZETA1= (L/10 . ) +T*VEL 

ZETA2= ( 9 . *L/10 . ) -T*VEL 

Kl= ( 3 . *MOD*INERT) / ( ZETA1**3 ) 

K2= ( 3 . *MOD*INERT) / ( ZETA2 **3-) 

STIFF= (K1*K2 ) / (K1+K2 ) 

DERY ( 1) = (P-STIFF*Y (2 ) ) /MASS 
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DERY ( 2 ) =Y ( 1 ) 

DEN0M=1.+ ( ZETA1**3/ZETA2**3) 

RETURN 

END 

ieicit’k'k’kie'kit'k'k'k'k’k'k'k’k'kie'k’k'k’k'kic'k'k’kic'kic'k’kit'k'k’kie'kic'kie'kic'k'k'k'kic'kicie'k'k 

****************************************************** 

SUBROUTINE OUTP (T, Y, DERY , IHLF , NDIM , PRMT ) 

DIMENSION Y ( 4 ) , DERY ( 4 ) 

COMMON P , L , VEL, MASS , INERT , K1 , K2 , STIFF , MOD , ZETA1 , ZETA2 , TF , DENOM 

DELTA2=Y ( 2 ) /DENOM 

DELTA1=Y ( 2 ) -DELTA2 

WRITE ( 15 , * ) T,Y(2) , DELTA 1 , DELTA2 

RETURN 

END 

****************************************************** 

****************************************************** 

C THE SUBROUTINE "RKGS" (RUNGE-KUTTA) IS INCLUDED 
C IN THE NEXT PROGRAM 

C RKG 

C RKG 

C RKG 

C SUBROUTINE RKGS RKG 

C RKG 

C PURPOSE RKG 

C TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL RKG 

C EQUATIONS WITH GIVEN INITIAL VALUES. RKG 

C RKG 

C USAGE RKG 

C CALL RKGS (PRMT, Y, DERY, NDIM, IHLF, FCT, OUTP, AUX) RKG 

C PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT. RKG 

C RKG 

C DESCRIPTION OF PARAMETERS RKG 

C PRMT - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER RKG 

C OR EQUAL TO 5, WHICH SPECIFIES THE PARAMETERS OF RKG 

C THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR RKG 

C COMMUNICATION BETWEEN OUTPUT SUBROUTINE (FURNISHED RKG 

C BY THE USER) AND SUBROUTINE RKGS. EXCEPT PRMT (5) RKG 

C THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE RKG 

C RKGS AND THEY ARE RKG 

C PRMT ( 1 ) - LOWER BOUND OF THE INTERVAL (INPUT), RKG 

C PRMT (2)- UPPER BOUND OF THE INTERVAL (INPUT), RKG 

100 : >TERVAL (INPUT), RKGS 220 

C PRMT (2)- UPPER BOUND OF THE INTERVAL 
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****** ******* *********************^********** 

V* ' ” 

C THIS PROGRAM SOLVES A SYSTEM OF 5 LINEAR * 

C DIFFERENTIAL EQUATIONS DESCRIBING THE * 

C DYNAMIC RESPONSE OF TWO COUPLED CANTILEVER * 

C BEAMS WITH A MOVING LOAD BETWEEN THEM. * 

C********************************************** 

c 

EXTERNAL FCT,OUTP 

DIMENSION FUNC(IO) , Y(8) , DERY ( 8 ) ,PRMT(5) ,AUX(8,8) ,A(5,5) 

REAL MOD , INERT , Ml , M2 , MB1 , MB2 ,31,82, OMEG , I PHI 1 , IPHI 2 
REAL I2PHI1,I2PHI2, LAMBDA, LENG 

COMMON IPHI1 , I PHI 2 , BETA, CONST , SI , S2 , Ml , M2 , MB1 , MB2 , VEL 
$ , LAMBDA, ZETA1 , ZETA2 , I2PHI1, I2PHI2 , LENG, PHI1 , PHI2 , DZZ 
WRITE (6 , *) 'ENTER APPLIED LOAD ( lbs ), AND SPEED OF MOVING LOAD' 
READ ( 5 , * ) P, VEL 

WRITE (6,*) 'ENTER BEAMLENGTH ( in) , AND BASE THICKNESS (in) ' 

READ (5,*) LENG, THICK 

WRITE ( 6 , * ) 'ENTER MASSES: Ml, M2, AND DENSITY OF BEAM ( lbs/ in3 ) ' 
READ ( 5 , * ) Ml, M2, DENS 
MOD=30 . E6 
LC=0 

INERT= (1 ./12 . ) *THICK**4 
DENS= (DENS/386 . 4 ) *THICK**2 


WRITE ( 6 , * ) 'APPLIED LOAD: ',P 

WRITE (6,*) 'VELOCITY : ' ,VEL 

WRITE ( 6 , * ) ' BEAM LENGTH : ' , LENG 

WRITE (6,*) 'THICKNESS : ', THICK 

WRITE (6,*) 'MASS/UNIT : ' , DENS 

C 


C**** ONE TERM OF THE NATURAL MODE SERIES IS USED ****** 

C 

OMEG= (1.875**2) * ( MOD* INERT/ ( DENS * LENG * * 4 ) ) * * 0 . 5 
BETA= ( ( (OMEG* *2 ) *DENS) / (MOD* INERT) ) **0.25 
WRITE ( 6 , * ) • OMEGA= ' , OMEG , ' BETA= ' , BETA 

C 

C**** EVALUATE THE CONSTANT TERMS IN MODE EQUATION ****** 

C 

X=BETA*LENG 
S1=SIN (X) -SINH(X) 

S2=COS (X) +COSH (X) 

WRITE ( 6 , *) 'S1=',S1,' S2=',S2 

C 

C**** EVALUATE CONSTANT OF INTEGRATION OF MODE EQUATION **** 
C**** INTEGRATE USING SIMPSON'S RULE - 100 INCREMENTS***** 

C 

XINC=LENG/100. 

ZETA=0 . 0 
XFACT=0 . 0 
DO 10 1=2,100,2 
DO 20 J=1 , 3 

ZETA=XINC* ( 1-2 ) +XINC* ( J-l) 

X=BETA*ZETA 

FUNC ( J) =S1**2* (SIN (X) **2+SINH (X) **2-2 . *SIN(X) *SINH(X) ) 
$+S2**2* (COS (X) **2+COSH (X) **2-2 . *COS (X) *COSH(X) ) 

$+2 . *S1*S2 * (SIN (X) *COS (X) -SIN(X) *COSH(X) -COS(X) *SINH(X) 
$+SINH (X) *COSH (X) ) 

20 CONTINUE 

XFACT=XFACT+ ( XINC/ 3 . ) * ( FUNC ( 1 ) +4 . *FUNC ( 2 ) +FUNC ( 3 ) ) 

10 CONTINUE 
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C0NST=-(1./(DENS*XFACT) ) **0.5 
WRITE ( 6 , * ) ' CONST= ' , CONST 

C 

C**** EVALUATE INTEGRAL OF MODE SHAPE: IPHI1=IPHI2**** 

C 

XINC=LENG/100. 

ZETA=0 . 0 
XFACT=0. 0 
DO 30 1=2,100,2 
DO 40 J=l, 3 

ZETA=XINC* (1-2 ) +XINC* ( J-l) 

X=ZETA*BETA 

FUNC ( J) = (SI* (SIN (X) -SINH (X) ) +S2* (COS (X) -COSH (X) ) ) 

40 CONTINUE 

XFACT=XFACT+ (XINC/3 . ) * (FUNC ( 1) +4 . *FUNC (2 ) +FUNC ( 3 ) ) 

30 CONTINUE 

IPHI1=DENS*XFACT*C0NST 

IPHI2=IPHI1 

WRITE (6,*) 1 IPHI1=IPHI2= 1 , IPHI1 

C 

C**** EVALUATE INTEGRAL OF (D2PHI/DZETA2 ) **2 : I2PHI **** 

C 

XINC=LENG/100. 

ZETA=0 . 0 
XFACT=0 . 0 
DO 50 1=2,100,2 
DO 60 J=l, 3 

ZETA=XINC* (1-2) +XINC* (J-l) 

X=Z ETA* BETA 

FUNC (J)= (Sl**2* (SIN (X) **2+SINH (X) **2 
$+2 . *SIN (X) *SINH(X) ) +S2**2* (COS (X) **2+COSH (X) **2 
$+2 . *COS (X) *COSH (X) ) +2 . *S1*S2* (SIN (X) *COS (X) +SIN (X) *COSH(X) 
$+COS(X) *SINH (X) +COSH (X) *SINH(X) ) ) 

60 CONTINUE 

XFACT=XFACT+ (XINC/3 . ) * (FUNC ( 1) +4 . *FUNC (2 ) +FUNC (3 ) ) 

50 CONTINUE 

1 2 PHI l=XFACT*MOD* INERT*CONST* * 2 * BETA* * 4 
I2PHI2=I2PHI1 

WRITE (6,*) ' I2PHI1=I2PHI2=' ,I2PHI1 

C 

MB1=DENS*LENG 
MB2=MB1 
ZETA1=LENG/10 . 

ZETA2=. 9*LENG 
X1=ZETA1*BETA 
X2=ZETA2*BETA 

PHI1=C0NST* (SI* (SIN (XI) -SINH (XI) ) +S2* (COS (XI) -COSH (XI) ) ) 
PHI2=CONST* (SI* (SIN (X2) -SINH (X2 ) )+S2* (COS (X2) -COSH (X2) ) ) 
WRITE ( 6 , * ) 1 PHI1= 1 , PHI1 , ' PHI2= ' , PHI2 

DZ=VEL 
DZZ=0 . 0 

D1PHI1=C0NST*BETA*DZ* (SI* (COS (XI) -COSH (XI) ) 

$-S2*(SIN(Xl) +SINH (XI) ) ) 

D1PHI2=C0NST*BETA* (-DZ) * (SI* (COS (X2 ) -COSH (X2 ) ) 

$-S2* (SIN (X2 ) +SINH (X2 ) ) ) 

C 

C**** DEFINE THE INITIAL VALUES OF THE DEPENDENT VARIABLES ** 

C 

Y (3 ) =P*PHI1/I2PHI1 
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Y (4 ) =P*PHI2/I2PHI2 
Y ( 1) =0 . 0 

Y (2 ) =Y (1) +PHI1*Y (3 ) +PHI2*Y (4 ) 

Y(7) =0 . 0 
Y (8 ) =0 . 0 
Y (5) =0 . 0 

Y(6) =Y (5) +Y (3) *D1PHI1+Y (4) *D1PHI2 
WRITE ( 6 , *) 'Y(1) = ',Y(1) 

WRITE ( 6 , * ) * Y ( 2 ) = ' , Y (2 ) 

WRITE ( 6 , *) 'Y(3)=',Y(3) 

WRITE ( 6 , * ) •Y(4)= , ,Y(4) 

DO 55 1=1,8 
DERY (I) =l./8 . 

55 CONTINUE ■ 

PRMT (1) =0 . 0 

PRMT ( 2 ) = . 8 *LENG/ ( VEL) 

PRMT ( 3 ) = ( PRMT ( 2 ) -PRMT (1))/100. 

PRMT (4) =0.0001 
NDIM=8 

CALL RKGS ( PRMT , Y , DERY , NDIM, IHLF , FCT , OUTP , AUX ) 
********************************************* 
********************************************* 

SUBROUTINE FCT(T, Y, DERY) 

DIMENSION Y (8) , DERY ( 8 ) ,A(5,5) ,B(5) 

REAL MOD , INERT , Ml , M2 , MB1 , MB2 , SI , S2 , OMEG , IPHI1 , 1 PHI 2 
REAL 1 2 PHI 1 , 1 2 PHI 2 , LAMBDA , LENG 

COMMON IPHI1, IPHI2 , BETA, CONST, SI , S2 ,M1,M2 ,MB1 , MB2 , VEL 
$ , LAMBDA, ZETA1 , ZETA2 , I2PHI1 , I2PHI2 , LENG, PHI1 , PHI2 , DZZ 
ZETA1= (LENG/10 . ) +VEL*T 
ZETA2= ( . 9*LENG) -VEL*T 
X1=ZETA1*BETA 
X2=ZETA2*BETA 
C 

C**** EVALUATE MODE SHAPE EQUATIONS AT EACH LOAD POSITION *** 
C 

PHI1=C0NST* (SI* (SIN (XI) -SINH (XI) ) +S2* (COS (XI) -COSH(Xl) ) ) 
PHI2=CONST* (SI* (SIN (X2 ) -SINH (X2) ) +S2* (COS (X2) — COSH (X2 ) ) ) 

C 

IF(LC.GT.IO) GO TO 71 
WRITE ( 6 , * ) ' **********TIME=' ,T 

C**** EVALUATE DERIVATIVES OF MODE SHAPE EQUATIONS **** 

C 

71 DZ=VEL 
DZZ=0 . 0 

D1PHI1=C0NST*BETA*DZ* (SI* (COS (XI) -COSH (XI) ) 

$-S2* (SIN (XI) +SINH(X1) ) ) 

D1PHI2=C0NST*BETA* (-DZ) * (SI* (COS (X2) -COSH(X2) ) 

$-S2* (SIN (X2) +SINH(X2) ) ) 

D2PHI1=C0NST*BETA* (SI* (-BETA*DZ**2*SIN (XI) +DZZ*COS (XI) 
$-BETA*DZ**2*SINH (XI) -DZZ*COSH (XI) ) -S2* (BETA*DZ**2*COS (XI) 
$+DZZ*SIN (XI) +BETA*DZ**2*COSH (XI) +DZZ*SINH (XI) ) ) 
D2PHI2=CONST*BETA* (SI* (-BETA*DZ**2*SIN(X2) +DZZ*COS (X2) 
$-BETA*DZ**2*SINH(X2)-DZZ*COSH(X2) )-S2*(BETA*DZ**2*COS(X2) 
$+DZZ*SIN(X2)+BETA*DZ**2*COSH(X2)+DZZ*SINH(X2) ) ) 

C 

C**** DEFINE MATRIX COEFFICIENTS A[5,5] **** 

C 

A ( 1 , 1 ) =M1+MB1 
A(2,l)=0.0 



A (3 , 1) =IPHI1 
A (4 , 1) =0 . 0 
A(5, 1)=-1. 

A(l,2)=0.0 
A ( 2 , 2 ) =M2+MB2 
A(3 / 2)=0.0 
A(4 , 2)=-IPHI2 
A ( 5 , 2 ) =1 . 0 
A ( 1 , 3 ) =IPHI1 
A(2,3)=0.0 
A(3,3)=1.0 
A(4,3)=0.0 
A (5 , 3 ) =-PHIl 
A ( 1 , 4 ) =0 . 0 
A (2 , 4 ) =-IPHI2 
A(3,4)=0.0 
A(4 / 4)=l. 0 
A(5,4)=-PHI2 
A(l,5)=1.0 
A(2 , 5) =-l . 0 
A ( 3 , 5) =PHI1 
A(4 , 5) =PHI2 
A(5,5)=0.0 
C 

C**** DEFINE RIGHT HAND SIDE B[5] ******* 

C 

B(l)=-P . 

B (2 ) =P 

B(3)=-I2PHI1*Y(3) 

B (4 ) =-I2PHI2*Y ( 4 ) 

B ( 5) =D2PHI1*Y (3 ) +D2PHI2*Y (4 ) +2 . *D1PHI1*Y (7 ) +2 . *D1PHI2*Y (8) 
IF(LC.GT.IO) GO TO 1 
DO 63 1=1,5 

63 WRITE (6, *) A (I , 1) ,A(I,2) ,A(I,3) ,A(I,4) ,A(I,5) ,B(I) 

1 LC=LC+1 

N=5 

CALL SIMQ (A, B, N, KS) 

IF(KS.EQ.l) WRITE (6, *) 'NO SOLUTION! !!!!! ' 

IF(LC.GT.IO) GO TO 64 

WRITE (6,*) 'BACK FROM SIMQ [B] ' 

DO 64 1=1,5 
WRITE ( 6 , * ) B(I) 

64 CONTINUE 
DERY (1) =Y ( 5 ) 

DERY ( 2 ) =Y (6) 

DERY (3) =Y (7) 

DERY (4) =Y ( 8 ) 

DERY (5) =B ( 1 ) 

DERY (6) =B ( 2 ) 

DERY (7) =B ( 3 ) 

DERY (8) =B(4) 

LAMBDA=B (5) 

IF (LC.GT. 10) GO TO 62 
DO 61 1=1,8 

61 WRITE (6 , *) 'Y(' ,1, •)=' ,Y(I) 

62 RETURN 
END 

************************************************** 

************************************************** 
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SUBROUTINE OUTP (T, Y, DERY , IHLF , NDIM , PRMT ) 

DIMENSION Y ( 8 ) , DERY ( 8 ) , PRMT ( 5 ) 

REAL MOD , INERT , Ml , M2 , MB1 , MB2 , SI , S2 , OMEG , IPHI1 , IPHI2 
REAL I 2 PHI 1,1 2 PHI 2, LAMBDA, LENG 

COMMON IPHI1 , I PH I 2 , BETA, CONST, SI , S2,M1,M2,MB1,MB2,VEL 
$, LAMBDA, ZETA1, ZETA2 , 12PHI1 , 12PHI2 , LENG, PHI 1 , PHI2 , DZZ 
ZETA1= (LENG/10 . ) +VEL*T 
ZETA2= ( . 9* LENG) -VEL*T 
X1=ZETA1*BETA 
X2=ZETA2*BETA 

PHIl=CONST* (SI* (SIN (XI) -SINH (XI) ) +S2* (COS (XI) -COSH (XI) ) ) 

PHI2=CONST* (SI* (SIN(X2) -SINH(X2) ) +S2* (COS (X2 ) -COSH (X2 ) ) ) 

WRITE (15,*) T, Y(2) -Y(l) , Y(3) *PHI1, Y (4) *PHI2 

WRITE (6,*) 1 ACCEL= 1 , DZZ 

RETURN 

END 

************************************************ 

************************************************ 

C 

c 

c 

C SUBROUTINE SIMQ 
C 

C PURPOSE 

C OBTAIN SOLUTION OF A SET OF SIMULTANEOUS LINEAR EQUATIONS, 

C AX=B 

C 

C USAGE 

C CALL SIMQ (A, B, N, KS) 

C 

C DESCRIPTION OF PARAMETERS 

C A - MATRIX OF COEFFICIENTS STORED COLUMNWISE. THESE ARE 

C DESTROYED IN THE COMPUTATION. THE SIZE OF MATRIX A IS 

C N BY N. 

C B - VECTOR OF ORIGINAL CONSTANTS (LENGTH N) . THESE ARE 

C REPLACED BY FINAL SOLUTION VALUES, VECTOR X. 

C N - NUMBER OF EQUATIONS AND VARIABLES. N MUST BE .GT. ONE. 

C KS - OUTPUT DIGIT 

C 0 FOR A NORMAL SOLUTION 

C 1 FOR A SINGULAR SET OF EQUATIONS 

C 

C REMARKS 

C MATRIX A MUST BE GENERAL. 

C IF MATRIX IS SINGULAR , SOLUTION VALUES ARE MEANINGLESS. 

C AN ALTERNATIVE SOLUTION MAY BE OBTAINED BY USING MATRIX 

C INVERSION (MINV) AND MATRIX PRODUCT (GMPRD) . 

C 

C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
C NONE 

C 

C METHOD 

C METHOD OF SOLUTION IS BY ELIMINATION USING LARGEST PIVOTAL 

C DIVISOR. EACH STAGE OF ELIMINATION CONSISTS OF INTERCHANGING 

C ROWS WHEN NECESSARY TO AVOID DIVISION BY ZERO OR SMALL 

C ELEMENTS . 

C THE FORWARD SOLUTION TO OBTAIN VARIABLE N IS DONE IN 

C N STAGES. THE BACK SOLUTION FOR THE OTHER VARIABLES IS 

C CALCULATED BY SUCCESSIVE SUBSTITUTIONS. FINAL SOLUTION 

C VALUES ARE DEVELOPED IN VECTOR B, WITH VARIABLE 1 IN B(l), 
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C VARIABLE 2 IN B(2), , VARIABLE N IN B(N) . 

C IF NO PIVOT CAN BE FOUND EXCEEDING A TOLERANCE OF 0.0, 

C THE MATRIX IS CONSIDERED SINGULAR AND KS IS SET TO 1. THIS 

C TOLERANCE CAN BE MODIFIED BY REPLACING THE FIRST STATEMENT. 

C 

c 

c 

SUBROUTINE SIMQ ( A, B, N , KS ) 

DIMENSION A( 1) ,B(1) 

FORWARD SOLUTION 

TOL=0 . 0 
KS=0 
JJ=-N 

DO 65 J=1 , N 
JY=J+1 
JJ=JJ+N+1 
BIGA=0 
IT=JJ-J 
DO 30 I=J,N 

SEARCH FOR MAXIMUM COEFFICIENT IN COLUMN 
IJ=IT+I 

IF ( ABS (BIGA) -ABS ( A ( IJ) ) ) 20,30,30 
20 BIGA=A(IJ) 

IMAX=I 
30 CONTINUE 

TEST FOR PIVOT LESS THAN TOLERANCE (SINGULAR MATRIX) 

IF (ABS (BIGA) -TOL) 35,35,40 
35 KS=1 
RETURN 

INTERCHANGE ROWS IF NECESSARY 

40 I1=J+N* ( J-2 ) 

IT=IMAX-J 
DO 50 K~J , N 
I1=I1+N 
I2=I1+IT 
SAVE=A(I1) 

A ( I 1 ) = A (12) 

A (12) =SAVE 

DIVIDE EQUATION BY LEADING COEFFICIENT 

50 A (II) =A(I1) /BIGA 
S AVE=B ( IMAX ) 

B (IMAX) =B ( J) 

B ( J) =SAVE/BIGA 

ELIMINATE NEXT VARIABLE 

IF(J-N) 55,70,55 
55 IQS=N*(J-1) 

DO 65 IX=JY , N 
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IXJ=IQS+IX 
IT=J-IX 
DO 60 JX=JY,N 
IXJX=N* (JX-1) +IX 
JJX=IXJX+IT 

60 A(IXJX) =A(IXJX) - (A(IXJ) *A(JJX) ) 
65 B(IX) =B(IX) - (B(J) *A(IXJ) ) 

C 

C BACK SOLUTION 

C 

70 NY=N— 1 
IT=N*N 

DO 80 J=1 , NY 
IA=IT-J 
IB=N-J 
IC=N 

DO 80 K=1 , J 

B (IB) =B (IB) -A (I A) *B(IC) 

IA=IA-N 
80 IC=IC-1 
RETURN 
END 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE RKGS 

PURPOSE 

TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL 

EQUATIONS WITH GIVEN INITIAL VALUES. 

USAGE 

CALL RKGS (PRMT, Y, DERY, NDIM, IHLF, FCT, OUTP, AUX) 

PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT. 

DESCRIPTION OF PARAMETERS 

PRMT - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER 
OR EQUAL TO 5, WHICH SPECIFIES THE PARAMETERS OF 
THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR 
COMMUNICATION BETWEEN OUTPUT SUBROUTINE (FURNISHED 
BY THE USER) AND SUBROUTINE RKGS. EXCEPT PRMT (5) 

THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE 
RKGS AND THEY ARE 

PRMT ( 1 ) - LOWER BOUND OF THE INTERVAL (INPUT), 

PRMT (2)- UPPER BOUND OF THE INTERVAL (INPUT), 

PRMT (3)- INITIAL INCREMENT OF THE INDEPENDENT VARIABLE 
(INPUT) , 

PRMT (4)- UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS 
GREATER THAN PRMT (4), INCREMENT GETS HALVED. 

IF INCREMENT IS LESS THAN PRMT (3) AND ABSOLUTE 
ERROR LESS THAN PRMT(4)/50, INCREMENT GETS DOUBLED. 
THE USER MAY CHANGE PRMT (4) BY MEANS OF HIS 
OUTPUT SUBROUTINE. 

PRMT ( 5 ) - NO INPUT PARAMETER. SUBROUTINE RKGS INITIALIZES 
PRMT (5) =0. IF THE USER WANTS TO TERMINATE 
SUBROUTINE RKGS AT ANY OUTPUT POINT, HE HAS TO 
CHANGE PRMT (5) TO NON-ZERO BY MEANS OF SUBROUTINE 
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c 
c 
c 
c 
c 
c 
c 

C Y 

C 

c 

C DERY 

C 
C 
C 

C NDIM 

C 

C IHLF 

C 
C 
C 
C 
C 
C 

C FCT 

C 
C 
C 
C 

C OUTP 

C 
C 
C 
C 
C 

C AUX 

C 
C 

C REMARKS 

C THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM, IF 

C (1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE 

C NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE 

C IHLF=11) , 

C (2) INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN 

C (ERROR MESSAGES IHLF=12 OR IHLF=13), 

C (3) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH, 

C (4) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO. 

C 

C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
C THE EXTERNAL SUBROUTINES FCT (X, Y , DERY) AND 

C OUTP(X,Y, DERY, IHLF, NDIM, PRMT) MUST BE FURNISHED BY THE USER 

C 

C METHOD 

C EVALUATION IS DONE BY MEANS OF FOURTH ORDER RUNGE-KUTTA 

C FORMULAE IN THE MODIFICATION DUE TO GILL. ACCURACY IS 

C TESTED COMPARING THE RESULTS OF THE PROCEDURE WITH SINGLE 

C AND DOUBLE INCREMENT. 

C SUBROUTINE RKGS AUTOMATICALLY ADJUSTS THE INCREMENT DURING 

C THE WHOLE COMPUTATION BY HALVING OR DOUBLING. IF MORE THAN 

C 10 BISECTIONS OF THE INCREMENT ARE NECESSARY TO GET 


OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE 
FEASIBLE IF ITS DIMENSION IS DEFINED GREATER 
THAN 5. HOWEVER SUBROUTINE RKGS DOES NOT REQUIRE 
AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL 
FOR HANDING RESULT VALUES TO THE MAIN PROGRAM 
(CALLING RKGS) WHICH ARE OBTAINED BY SPECIAL 
MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP. 

- INPUT VECTOR OF INITIAL VALUES. (DESTROYED) 
LATERON Y IS THE RESULTING VECTOR OF DEPENDENT 
VARIABLES COMPUTED AT INTERMEDIATE POINTS X. 

- INPUT VECTOR OF ERROR WEIGHTS. (DESTROYED) 

THE SUM OF ITS COMPONENTS MUST BE EQUAL TO 1. 
LATERON DERY IS THE VECTOR OF DERIVATIVES, WHICH 
BELONG TO FUNCTION VALUES Y AT A POINT X. 

-AN INPUT VALUE, WHICH SPECIFIES THE NUMBER OF 
EQUATIONS IN THE SYSTEM. 

- AN OUTPUT VALUE, WHICH SPECIFIES THE NUMBER OF 
BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS 
GREATER THAN 10, SUBROUTINE RKGS RETURNS WITH 
ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM. ERROR 
MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE 

PRMT ( 3 ) =0 OR IN CASE SIGN ( PRMT ( 3 ) ) . NE . SIGN ( PRMT ( 2 ) 
PRMT ( 1 ) ) RESPECTIVELY. 

- THE NAME OF AN EXTERNAL SUBROUTINE USED. THIS 
SUBROUTINE COMPUTES THE RIGHT HAND SIDES DERY OF 
THE SYSTEM TO GIVEN VALUES X AND Y. ITS PARAMETER 
LIST MUST BE X, Y , DERY . SUBROUTINE FCT SHOULD 

NOT DESTROY X AND Y. 

- THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED. 

ITS PARAMETER LIST MUST BE X,Y, DERY, IHLF, NDIM, PRMT 
NONE OF THESE PARAMETERS (EXCEPT, IF NECESSARY, 
PRMT (4) , PRMT (5) , . . . ) SHOULD BE CHANGED BY 
SUBROUTINE OUTP. IF PRMT (5) IS CHANGED TO NON-ZERO 
SUBROUTINE RKGS IS TERMINATED. 

- AN AUXILIARY STORAGE ARRAY WITH 8 ROWS AND NDIM 
COLUMNS . 
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SATISFACTORY ACCURACY, THE SUBROUTINE RETURNS WITH 
ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM. 

TO GET FULL FLEXIBILITY IN OUTPUT, AN OUTPUT SUBROUTINE 
MUST BE FURNISHED BY THE USER. 

FOR REFERENCE, SEE 

RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL COMPUTERS, 
WILEY, NEW YORK/LONDON, 1960, PP. 110-120. 


SUBROUTINE RKGS ( PRMT , Y , DERY , NDIM , IHLF , FCT , OUTP , AUX ) 


DIMENSION Y ( 1) , DERY ( 1 ) ,AUX(8,1) ,A(4) ,B(4) ,C(4) ,PRMT(1) 
DO 1 1=1, NDIM 

1 AUX (8 , I) = . 06666667*DERY (I) 

X=PRMT ( 1 ) 

XEND=PRMT ( 2 ) 

H=PRMT ( 3 ) 

PRMT (5) =0. 

CALL FCT (X,Y, DERY) 

ERROR TEST 

IF (H* (XEND-X) ) 38 , 37 , 2 

PREPARATIONS FOR RUNGE-KUTTA METHOD 

2 A(l) =. 5 
A(2)=. 2928932 
A(3)=l. 707107 
A(4)=. 1666667 
B ( 1) =2 . 

B (2 ) =1 . 

B (3 ) =1 . 

B (4 ) =2 . 

C (1) =• 5 
C(2)=. 2928932 
C(3)=l. 707107 
C (4 ) = . 5 

PREPARATIONS OF FIRST RUNGE-KUTTA STEP 
DO 3 1=1, NDIM 
AUX (1,1) =Y ( I ) 

AUX (2,1) =DERY ( I ) 

AUX (3, I) =0. 

3 AUX (6, I) =0. 

IREC=0 
H=H+H 
IHLF=— 1 
ISTEP=0 
IEND=0 


START OF A RUNGE-KUTTA STEP 

4 IF ( (X+H-XEND) *H) 7,6,5 

5 H=XEND-X 

6 IEND=1 

RECORDING OF INITIAL VALUES OF THIS STEP 

7 CALL OUTP (X,Y, DERY, IREC, NDIM, PRMT) 
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IF ( PRMT (5))40,8,40 

8 ITEST=0 

9 ISTEP=ISTEP+1 


START OF INNERMOST RUNGE-KUTTA LOOP 
J=1 

10 AJ=A( J) 

BJ=B(J) 

CJ=C(J) 

DO 11 1=1 , NDIM 
R1=H*DERY ( I ) 

R2=AJ * ( Rl-BJ * AUX (6,1) ) 

Y(I)=Y(I) +R2 
R2=R2+R2+R2 

11 AUX (6 , I ) =AUX (6,1) +R2-CJ*R1 
IF ( J-4 ) 12,15,15 

12 J=J+1 

IF ( J-3 ) 13,14,13 

13 X=X+.5*H 

14 CALL FCT (X , Y , DERY) 

GOTO 10 

END OF INNERMOST RUNGE-KUTTA LOOP 


TEST OF ACCURACY 

15 IF (ITEST) 16,16,20 

IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY 

16 DO 17 1=1, NDIM 

17 AUX (4,1) =Y ( I ) 

ITEST=1 

ISTEP=ISTEP+ISTEP-2 

18 IHLF=IHLF+1 
X=X-H 
H=.5*H 

DO 19 1=1, NDIM 
Y ( I ) =AUX (1,1) 

DERY ( I ) =AUX (2,1) 

19 AUX (6,1) =AUX (3,1) 

GOTO 9 

IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE 

20 IMOD=ISTEP/2 

IF ( ISTEP-IMOD-IMOD) 21,23,21 

21 CALL FCT (X , Y , DERY) 

DO 22 1=1, NDIM 
AUX (5,1) =Y ( I ) 

22 AUX (7,1) =DERY ( I ) 

GOTO 9 

COMPUTATION OF TEST VALUE DELT 

23 DELT=0 . 

DO 24 1=1, NDIM 

24 DELT=DELT+AUX (8,1) *ABS (AUX (4,1) -Y (I) ) 

IF(DELT-PRMT(4) )28,28,25 


ERROR IS TOO GREAT 
25 IF (IHLF-10) 26,36,36 
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nnn non no 


26 DO 27 1=1 , NDIM 

27 AUX (4,1) =AUX (5,1) 

ISTEP=ISTEP+ISTEP-4 

X=X-H 

IEND=0 
GOTO 18 

RESULT VALUES ARE GOOD 

28 CALL FCT ( X , Y , DERY ) 

DO 29 1=1, NDIM 
AUX (1,1) =Y ( I ) 

AUX (2,1) =DERY ( I ) 

AUX (3,1) =AUX (6,1) 

Y ( I ) =AUX (5,1) 

29 DERY ( I ) =AUX (7,1) 

CALL OUTP ( X-H , Y , DERY , IHLF , NDIM , PRMT ) 

IF (PRMT(5) )40,30,40 

30 DO 31 1=1, NDIM 
Y ( I ) =AUX (1,1) 

31 DERY ( I ) =AUX (2,1) 

IREC=IHLF 

IF (IEND) 32,32,39 

INCREMENT GETS DOUBLED 

32 IHLF=IHLF-1 

Tcmu t> — rcmpD / o 

XUliil — XU 

H=H+H 

TO STOP AT INPUT DELT: IF (IHLF) 4,33,33 
IF (IHLF) 4,33,33 

33 IMOD=ISTEP/2 

IF ( ISTEP-IMOD-IMOD) 4,34,4 

34 IF(DELT-. 02*PRMT (4 )) 35,35,4 

35 IHLF=IHLF-1 
ISTEP=ISTEP/2 
H=H+H 

GOTO 4 


RETURNS TO CALLING PROGRAM 

36 IHLF=11 

CALL FCT (X,Y, DERY) 

GOTO 39 

37 IHLF=12 
GOTO 39 

38 IHLF=13 

39 CALL OUTP (X,Y, DERY, IHLF, NDIM, PRMT) 

40 RETURN 
END 
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